{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 2D Optimal transport for different metrics\n",
    "\n",
    "\n",
    "2D OT on empirical distributio  with different gound metric.\n",
    "\n",
    "Stole the figure idea from Fig. 1 and 2 in\n",
    "https://arxiv.org/pdf/1706.07650.pdf\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Remi Flamary <remi.flamary@unice.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "import ot.plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset 1 : uniform sampling\n",
    "----------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAADSCAYAAAAWl/SpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFV1JREFUeJzt3X+0bGV93/H3J4C/QAvILaKCV4WY0FTRdUCyFm0gagLGVUxqKLS48FdJqLba0twSf8AVY7Ny2xKbttFiRQgY9EaN0gS7IHgMmhrgoKAoISIBgVwuRwgCapEf3/6x95Vh7j33zD17Zp8z57xfa82amb337OeZ58yZ7zx7f5/9pKqQJGml+4nlroAkSaMwYEmSpoIBS5I0FQxYkqSpYMCSJE0FA5YkaSoYsLTmJXlDki8tdz3GKUklObh9/KEk7xnTfg9K8mCS3drnX0jylnHsu93f55KcMq79aXUxYIkkRyX5v0m+l+TeJH+R5PDlrtdKkGR9++W/+zLW4dYkr1zq66vq16vqfeMop6q+U1V7VdWjS63PQHkbk1w0tP/jquqCrvvW6rRs/4RaGZI8A/gT4DRgM/Ak4B8BD02grN2r6pFx73clW03veTW9F00ne1j6SYCquriqHq2qH1bVZVX1NYAkP5Hk3UluS3J3kj9I8vfadUcnuWNwZ4O/0ttf0J9MclGS+4E3JNktyTuTfDvJA0muTXJgu/1PJbm87eXdlOSEhSqd5I1Jbmz3cUuSXxtYd3SSO5Kc3tZ5S5I3Dqx/ZpJLktyf5GrghTtpnyvb+/vaQ2E/m+SFST6f5J4k303ysSR7D7XBf0jyNeD7SXZP8rIkX23r+0dJPpHktwZe85ok1yW5r+3tvrhdfiFwEPC/2/I3LNAev9G+z79N8qahdedvKyvJfkn+pC3n3iRfbP/G25Uz0Lt8c5LvAJ9foMf5wiRXt+352ST7Dv4dhupya5JXJjkWeCfwz9ryrm/X//gQ4yKfvW31OCXJd9q/w7sGyjkiyVxbp61JztnJ31jToqq8reEb8AzgHuAC4Dhgn6H1bwJuBl4A7AV8GriwXXc0cMfQ9rcCr2wfbwQeBl5L8+PoqcBvAF8HXgQEeAnwTGBP4HbgjTQ9/5cC3wUOXaDev0QTaAL8HPAD4GUD9XoEOBvYA3h1u36fdv3HaXqTewI/A9wJfGmBctYDBew+sOxg4FXAk4F1NEHtA0NtcB1wYPuenwTcBry9rc+vAD8Cfqvd/qXA3cDLgd2AU9p9PHm4TReo47HA1va97An8YVvng9v15w+U9dvAh9p67EHTm86Oyhl473/Q7vepw+0BfKFtv21lfwq4aBc+HxcNrf8C8JYRPnvb6vHhtl4voTkq8NPt+i8Dr28f7wUcudz/a9663+xhrXFVdT9wFI//88+3vY/9203+BXBOVd1SVQ8CvwmcmNHP6Xy5qj5TVY9V1Q+BtwDvrqqbqnF9Vd0DvAa4tao+WlWPVNVXab78fnWBev9pVX273cefA5fRfPlu8zBwdlU9XFWXAg8CL0qTLPBPgTOr6vtVdQNNsB5ZVd1cVZdX1UNVNQ+cQxM0B/1eVd3evucjaYLw77X1+TRw9cC2pwL/s6quqqaXewHNl++RI1bpBOCjVXVDVX2fJhAs5GHgAOB5bV2+WFWLXVB0Y9tWP1xg/YUDZb8HOKFt565G+ey9t5qjAtcD19MELmje58FJ9quqB6vqL8dQHy0zA5aoqhur6g1V9VyaX8rPBj7Qrn42Te9gm9tovnz3ZzS3Dz0/EPj2DrZ7HvDy9lDVfUnuo/nCetaOdprkuCR/2R7Wuo+mF7XfwCb31BPPt/yA5pf2urb+g/UafH+LSrJ/ko8nubM91HnRUNkM7f/ZwJ1DgWFw/fOA04fe+4Ht60bxbEZ/P/+JptdyWXso9YwR9j/8N9zZ+ttoem7D7bEUo3z27hp4vO1vDPBmmsPdf5XkmiSvGUN9tMwMWHqCqvormkNIP9Mu+luaL9RtDqI53LYV+D7wtG0r2l/V64Z3OfT8dnZ8zuh24M+rau+B215VddrwhkmeTNP7+s/A/lW1N3ApzeHBxcy39T9w6D0tZEe9j//YLv+HVfUM4OQdlD34ui3Ac5IMbjNY/u3A+4fe+9Oq6uKd1GHQFkZ8P1X1QFWdXlUvAP4J8O+SvGKRchYrf7jsh2kO5y72+Vhsvzv77O1UVX2rqk4C/j7wO8Ank+y52Ou0shmw1rg0iQ6nJ3lu+/xA4CRg2yGUi4F/m+T5Sfai+bL+RNt7+WvgKUl+KckewLtpzuvszP8C3pfkkDRenOSZNJmKP5nk9Un2aG+HJ/npHezjSW0588AjSY4DfmGU91tNOvangY1JnpbkUJpzRguZBx6jOY+yzdNpDjF+L8lzaM7L7cyXgUeBt7UJGMcDRwys/zDw60le3rbJnm2bPr1dv3Wo/GGbaRJaDk3yNOCshTZMk9xxcBs8v9fW67ERy1nIyQNlnw18sm3nxT4fW4H1SRb6HtrZZ2+nkpycZF1VPQbc1y5+bGev0cpnwNIDNCf7r0ryfZpAdQNwerv+POBCmsSCvwH+H/CvAarqe8C/oglCd9L8on5CVtgOnEPzBXsZcD/wEeCpVfUATdA5keaX9V00v4y3C4Dttv+m3c/fAf8cuGQX3vPbaA4d3UXTm/zoQhtW1Q+A9wN/0R6uOxJ4L/Aymi/8P6UJgAuqqh/RJFq8mebL82SaAP1Qu34O+JfAf2/fz83AGwZ28dvAu9vy//0O9v85mkO4n29f+/mdVOcQ4M9oAu6Xgd+vqtlRytmJC2na8S7gKTR/m1E+H3/U3t+T5Cs72O+Cn70RHAt8I8mDwH8FTtzJOThNiW3ZQZJ6lOQq4ENVtWCwlPRE9rCkHiT5uSTPag8JngK8GPg/y10vaZp4pQupHy/i8bFftwCvq6oty1slabp4SFCSNBU8JChJmgoGLEnSVOj1HNZ+++1X69ev77NISdIKd+211363qoYvOrCdXgPW+vXrmZub67NISdIKl2Sky6N5SFCSNBUMWJKkqbBowEpyYJLZJN9M8o0kb2+X75tmsr1vtff7TL66WpE2bYLZ2Scum51tlkvSmIzSw3oEOL2qDqWZn+et7QVDzwCuqKpDgCva51qLDj8cTjjh8aA1O9s8P/zw5a2XpFVl0YBVVVuq6ivt4weAG4HnAMfz+MR3F9DMKqu16JhjYPPmJkideWZzv3lzs1ySxmSXzmElWU8znfdVNPMQbbu0zF0sMKFfklOTzCWZm5+f71BVrWjHHAOnnQbve19zb7CSNGYjB6x2PppPAe9op1X/sXYm1R1e46mqzq2qmaqaWbdu0TR7TavZWfjgB+E972nuh89pSVJHIwWsdvK1TwEfq6ptc/9sTXJAu/4A4O7JVFEr3rZzVps3w9lnP3540KAlaYxGyRIMzSR7N1bVOQOrLuHxmVpPAT47/uppKlxzzRPPWW07p3XNNctbL0mryqJXa09yFPBF4Os8PsX0O2nOY20GDgJuA06oqnt3tq+ZmZnySheSpEFJrq2qmcW2W/TSTFX1JSALrH7FrlZMq9CmTU0K+2Cixexs08PasGH56iVpVfFKF+rOcViSeuCMw+pucBzWaac1WYKOw5I0ZvawNB6Ow5I0YQYsjYfjsCRNmAFL3TkOS1IPDFjqznFYknpgwJIkTQUDlrozrV1SD0xrV3emtUvqgT0sjYdp7ZImzICl8TCtXdKEGbDUnWntknpgwFJ3prVL6oEBS5I0FQxY6s60dkk9MK1d3ZnWLqkH9rA0Hqa1S5owA5bGw7R2SRNmwFJ3prVL6oEBS92Z1i6pB6mq3gqbmZmpubm53sqTJK18Sa6tqpnFtrOHpe42bdr+8N/sbLNcksbEgKXuHIclqQeOw1J3jsOS1AN7WBoPx2FJmjADlsbDcViSJsyApe4chyWpBwYsdec4LEk9WDRgJTkvyd1JbhhYtjHJnUmua2+vnmw1JUlr3Sg9rPOBY3ew/Her6rD2dul4q6WpYlq7pB4sGrCq6krg3h7qomk1mNZ+5pmPn88yU1DSGHU5h/W2JF9rDxnus9BGSU5NMpdkbn5+vkNxWtFMa5c0YUsNWB8EXggcBmwB/stCG1bVuVU1U1Uz69atW2JxWvFMa5c0YUsKWFW1taoerarHgA8DR4y3WpoqprVL6sGSAlaSAwae/jJww0Lbag0wrV1SDxadXiTJxcDRwH7AVuCs9vlhQAG3Ar9WVVsWK8zpRSRJw0adXmTRi99W1Uk7WPyRJdVKq9OmTU0K+2Cixexs08PasGH56iVpVfFKF+rOcViSeuD0IurO6UUk9cAelsbDcViSJsyApfFwHJakCTNgqTvHYUnqgQFL3TkOS1IPDFiSpKlgwFJ3prVL6oFp7erOtHZJPbCHpfEwrV3ShBmwNB6mtUuaMAOWujOtXVIPDFjqzrR2ST0wYEmSpoIBS92Z1i6pB6a1qzvT2iX1wB6WxsO0dkkTZsDSeJjWLmnCDFjqzrR2ST0wYKk709ol9SBV1VthMzMzNTc311t5kqSVL8m1VTWz2Hb2sNTdpk3bH/6bnW2WS9KYGLDUneOwJPXAcVjqznFYknpgD0vj4TgsSRNmwNJ4OA5L0oQZsNSd47Ak9cCApe4chyWpB4sGrCTnJbk7yQ0Dy/ZNcnmSb7X3+0y2mpKktW6UHtb5wLFDy84ArqiqQ4Ar2udaq0xrl9SDRQNWVV0J3Du0+HjggvbxBcBrx1wvTZPBtPYzz3z8fJaZgpLGaKnnsPavqi3t47uA/RfaMMmpSeaSzM3Pzy+xOK14prVLmrDOSRfVXIxwwQsSVtW5VTVTVTPr1q3rWpxWKtPaJU3YUgPW1iQHALT3d4+vSpo6prVL6sFSA9YlwCnt41OAz46nOppKprVL6sGi04skuRg4GtgP2AqcBXwG2AwcBNwGnFBVw4kZ23F6EUnSsFGnF1n04rdVddICq16xy7XS6rRpU5PCPphoMTvb9LA2bFi+eklaVbzShbpzHJakHji9iLpzehFJPbCHpfFwHJakCTNgaTwchyVpwgxY6s5xWJJ6YMBSd47DktQDA5YkaSoYsNSdae2SemBau7ozrV1SD+xhaTxMa5c0YQYsjYdp7ZImzICl7kxrl9QDA5a6M61dUg8MWJKkqWDAUnemtUvqgWnt6s60dkk9sIel8TCtXdKEGbA0Hqa1S5owA5a6M61dUg8MWOrOtHZJPUhV9VbYzMxMzc3N9VaeJGnlS3JtVc0stp09LHW3adP2h/9mZ5vlkjQmBix15zgsST1wHJa6cxyWpB7Yw9J4OA5L0oQZsDQejsOSNGEGLHXnOCxJPTBgqTvHYUnqQaekiyS3Ag8AjwKPjJJHL0nSUoyjh3VMVR1msFrDTGuX1APT2tWdae2SetC1h1XAZUmuTXLqjjZIcmqSuSRz8/PzHYvTimVau6QJ6xqwjqqqlwHHAW9N8o+HN6iqc6tqpqpm1q1b17E4rVimtUuasE4Bq6rubO/vBv4YOGIcldKUMa1dUg+WHLCS7Jnk6dseA78A3DCuimmKmNYuqQdLnl4kyQtoelXQJG/8YVW9f2evcXoRSdKwUacXWXKWYFXdArxkqa/XKrJpU5PCPphoMTvb9LA2bFi+eklaVbzShbpzHJakHjgOS905DktSD+xhaTwchyVpwgxYGg/HYUmaMAOWunMclqQeGLDUneOwJPXAgCVJmgoGLHVnWrukHpjWru5Ma5fUA3tYGg/T2iVNmAFL42Fau6QJM2CpO9PaJfXAgKXuTGuX1AMDliRpKhiw1J1p7ZJ6YFq7ujOtXVIP7GFpPExrlzRhBiyNh2ntkibMgKXuTGuX1AMDlrozrV1SD1JVvRU2MzNTc3NzvZUnSVr5klxbVTOLbWcPS91t2rT94b/Z2Wa5JI2JAUvdOQ5LUg8ch6XuHIclqQf2sDQejsOSNGEGLI2H47AkTZgBS905DktSDwxY6s5xWJJ60ClgJTk2yU1Jbk5yxrgqpSmzYcP256yOOaZZPmTjxoV3s9R1y/Va67S8+9Xas+SBw0l2A/4aeBVwB3ANcFJVfXOh1zhwWAks9JFb6rrleq11Wv46aXXoY+DwEcDNVXVLVf0I+DhwfIf9SZK0oC4B6znA7QPP72iXPUGSU5PMJZmbn5/vUJym1caNzS/lpHm+7fHGjUtf12W/1mm666S1q8shwdcBx1bVW9rnrwdeXlVvW+g1HhLUSjysZJ2mt05aHfo4JHgncODA8+e2yyRJGrsuAesa4JAkz0/yJOBE4JLxVEur1VlnjX/dcr3WOi1/nbS2dJpeJMmrgQ8AuwHnVdX7d7a9hwQlScNGPSTY6eK3VXUpcGmXfUiSNAqvdCFJmgoGLEnSVOh0DmuXC0vmgdt6K3C89gO+u9yVmAK202hsp9HYTqOZ9nZ6XlWtW2yjXgPWNEsyN8pJwbXOdhqN7TQa22k0a6WdPCQoSZoKBixJ0lQwYI3u3OWuwJSwnUZjO43GdhrNmmgnz2FJkqaCPSxJ0lQwYC0iya8m+UaSx5LMDK37zXa25ZuS/OJy1XGlcAbqHUtyXpK7k9wwsGzfJJcn+VZ7v89y1nElSHJgktkk32z/597eLretBiR5SpKrk1zfttN72+XPT3JV+//3ifYar6uKAWtxNwC/Alw5uDDJoTQX/P0HwLHA77ezMK9J7Xv/H8BxwKHASW0bCc6n+YwMOgO4oqoOAa5on691jwCnV9WhwJHAW9vPkG31RA8BP19VLwEOA45NciTwO8DvVtXBwN8Bb17GOk6EAWsRVXVjVd20g1XHAx+vqoeq6m+Am2lmYV6rnIF6AVV1JXDv0OLjgQvaxxcAr+21UitQVW2pqq+0jx8AbqSZFNa2GlCNB9une7S3An4e+GS7fFW2kwFr6UaacXkNsT12zf5VtaV9fBew/3JWZqVJsh54KXAVttV2kuyW5DrgbuBy4NvAfVX1SLvJqvz/63S19tUiyZ8Bz9rBqndV1Wf7ro/WlqqqJKbrtpLsBXwKeEdV3Z/kx+tsq0ZVPQoclmRv4I+Bn1rmKvXCgAVU1SuX8DJnXH4i22PXbE1yQFVtSXIAzS/lNS/JHjTB6mNV9el2sW21gKq6L8ks8LPA3kl2b3tZq/L/z0OCS3cJcGKSJyd5PnAIcPUy12k5OQP1rrkEOKV9fAqw5nvyabpSHwFurKpzBlbZVgOSrGt7ViR5KvAqmvN9s8Dr2s1WZTs5cHgRSX4Z+G/AOuA+4Lqq+sV23buAN9FkN72jqj63bBVdAXZ1Buq1IsnFwNE0V9TeCpwFfAbYDBxEM4PBCVU1nJixpiQ5Cvgi8HXgsXbxO2nOY9lWrSQvpkmq2I2m07G5qs5O8gKaZKd9ga8CJ1fVQ8tX0/EzYEmSpoKHBCVJU8GAJUmaCgYsSdJUMGBJkqaCAUuSNBUMWJKkqWDAkiRNBQOWJGkq/H+2HAKLMVcUowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAACrCAYAAACZks5mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXuULVV95z/f8+7u+wZEnmJ8JIIzshSfgwbfQuJgJi8xMWBM0Kw4jiu+iGs5ksQYklFxMhgNKrmoESWTYPAtShjGjA9gBhXxGbwIeAF5KSDK7e49f1S1nHu6ftWnzjl9+hT3+1mrV5+zd+3a9fjV/lXV+f5+WykljDHGGFMvGhu9AcYYY4ypjh24McYYU0PswI0xxpgaYgdujDHG1BA7cGOMMaaG2IEbY4wxNcQOvA9Jx0m6vu/71yQdN8yyxqwnkpKkh06xvyPyPlv5909IOnmYZY2pgqTTJb0//3y4pLskNdda1tTYgUvaJeme/GSv/J01yT5SSkellC6Z5DpnDUmnSPrcRm/HNJB0rKT/I+mHkm6T9K+SHrvR2zUuki6R9JOBa+Ejk+wjpXR8SuncSa5z1qjbTfmk7Tm3o98rqV+5Ubtr4O83R+1zkJTS91JKm1JKS5Na5ywiaaekN467nrrfMT83pfSZjd4IM/tI2gJ8FPgD4HygAzwZ+OkGbEtzHQaol6WU3j3hdZoZZZL2LEmAKjTZllJarNqPmTy1fQIvY/A1S8HrwB2S/k7S9yXdLunDwXp2SXpG/nkuv2u6XdLVwGMHlj1Y0j9K+oGk70p6eV/d4yR9XtIdknZLOktSp68+SXqppG/ny7w9v6iKtqkp6XWS/k3SnZKukHRYXvckSZfld+SXSXpSX7tTJF2Tt/mupN+S9AjgncAT8zvpO0Y43HXh4QAppfNSSksppXtSSp9OKX0FfnZc3yzplvw4/eGAzfzMFvLvgzb2D5JuzI/9pZKO6qvbKekdkj4u6W7gqZK6eX/fk3STpHdKmutr8+rcVr4v6XdH3emiNyzqex2f2/VbJF2bb/vn+rejr83Pns4GjxXwSwPLbpX0nnz7b5D0RuWvRCU9RNLFkm7N2/+9pG19bXdJepWkr+Tb8yFJvZL9+31JX8/t+mpJj87LH5Fv8x3Kfgr7j31tTsiXvTPfvldJWgA+ARys+54sDx7hkE+Lce35Ekl/LulfgR8D7yO7AThLI77N1MAT/KDtSTpK0kXK3hbcJOl1BesYHKsfLOl/5efqImD/geWfoOwtxB2Svqy+nzwlvajPNq6R9JK+uuMkXS/plZJuzm31RSX7FvqM3Aa/k+/XhSt2o4wz8/X/SNJXJT1S0qnAbwGv0bhvy1JKtfwDdgHPCOpOB97f9/0IIAGt/PvHgA8B24E28It5+XHA9UV9AGcA/xvYARwGXLWyLNmN0BXAfyW7E/454Brg2Xn9Y4AnkL3xOAL4OvCKvn4S2d30NuBw4AfAc4J9ezXwVeDnye6aHwXsl2/X7cAL835Oyr/vBywAPwJ+Pl/HQcBR+edTgM9t9Pmcgr1sAW4FzgWOB7YP1L8U+EZ+bncA/zJgM3vZW4GN/S6wGegCbwOu7KvbCfwQ+A+5rfSAM4EL8742Ax8B/iJf/jnATcAj83P3gXxbHhrs2yXA7wV1q85v/7qAt+ftDwGawJPyfThiYP9/1scQx+oC4G/zbX8A8CXgJXndQ4Fn5n0cAFwKvG3gmvsScHC+7q8DLw327deBG8huppWv+0Fk1/R3gNeRXY9PA+7ss//dwJPzz9uBRxdd/7P8x/j2fAnwPeAosvGiXWZHeZu9bGItO+y3PTIb3w28ksz+NwOPH7yWCuzu88Bbc3t5Sn4eV5Y9JD8GJ5BdV8/Mvx+Q1/8S8JDcNn6R7Eal/1wvAn+a7/sJef32YN8in/E04Bbg0fk2/g/g0rzu2WR+YVu+DY8ADuobE944th1stCGOYcC7gLuAO/r+fn/QIAaNgsx5LRedKMod+DX0OVXgVO5z4I8Hvjewrj8G/i7Y9lcAF/R9T8Cxfd/PB04L2n4TOLGg/IXAlwbKPk92ES3kx+dXgbmBZU5hH3Dg+b4+Ir9wrs8v3guBA/O6i+lzFMCzqODAB/rZlrfdmn/fCby3r17A3cBD+sqeCHw3/3wOcEZf3cNZ24H/eOBa+LPo/K6si2zQuwd4VME6f3bN9PWx4sDDYwUcSPYad66v/iTgX4Jtfx7w//q+7wJ+u+/7XwHvDNp+CvgvBeVPBm4EGn1l5wGn55+/B7wE2DLQ7jhq4sAnYM+XAH9aYEfDOPA7Bv4eUdSevR34Sf3neWC9p1PgwMkeZhaBhb5lP9C37GuB9xXYxMlBPx9esZf8XN9D380IcDPwhIJ2ZT7jPcBf9X3fBOzJ9+NpwLfIHt4aA+12MgEHXvdX6M9LKW3r+3vXEG0OA25LKd1esa+Dgev6vl/b9/lBZK/e7lj5I7v7PxBA0sMlfVTZK9YfAW9i4FUQ2YCzwo/JDCHa/n8Ltu/agbJrgUNSSncDv0l2V75b0sck/UK0o/dXUkpfTymdklI6lOzp9mCyp2UoP7+l5K8rz1D2s8aPyJwQ7H2O+9d9ADAPXNFnL5/My0fdlpcPXAuvH6LN/mRPQ0X2VMZa10KbzM5W9u1vyZ7EkXSgpA/mr65/BLyf9bkWrkspLQ9s4yH5518le+K6Nn89+8Rg/TPNBOz5uoKyYdh/wNa+PkSb6FyVcTBwez5+rTBoa78+MO4eS+ZwkXS8pC/kr7bvIDvn/bZ2a9r7t/zI1sp8xl7jbkrpLrK3AIeklC4GziJ7y3WzpLOVaRcmRt0deMTdZAPkCg/s+3wdsEN9v7sNyW6yE7nC4QPr/O6AUW9OKZ2Q17+D7HXWw1JKW8icexXRSD/Xkb0WGuT7ZAbdz+FkrxhJKX0qpfRMMuP+BrBys5NG3I5ak1L6Btld8CPzorLzC+U29QLgROAZwFayu2/Y+xz3H+dbyO7+j+qzl60ppZXBY61tqcJe2y2pf7tvAX5CsT2Vsda18FP2HuS3pJRWNAFvIjsW/y6/Fn6b9bkWDpPUP771XwuXpZROJLup+DDZGy+o8bUwgj3D6v0dd//XGnd/ruL6dgPbc33CCoO29r6BcXchpXSGpC7wj8Cbyd5KbAM+zmi2VuYz9hp3823dj/ts7a9TSo8BjiR7k/bqfNGJ2Nr91YFfCTxFWUzhVrLX2QCklHaTiVX+RtJ2SW1JTxlinecDf5y3ORT4z311XwLulPRaZaKgZi5WWBG6bSb7Dfqu/Mn3D8bYt3cDfybpYblI4t9L2o/MOB8u6QWSWspCO44EPpo/9ZyYG9dPyX56WHk6uQk4VH2iuvsjkn4hF6wcmn8/jOy13hfyRc4HXi7pUEnbgdMGVnEl8PzcXo4Bfq2vbjPZcb2VbAB7U9m25E+G7wLOlLTyZHqIpGf3bcspko6UNA+8YbS9BuDLwFGSjlYmBjt9YDvOAd6qTITZlPTEfPArIzxW+fX1aeAtkrZIaigTrv1ivshmMvv7oaRDuG9AG4V3A6+S9Jj8WniopAcBXyR7mnpNfr6OA54LfFBSR5mAc2tKaQ/Zddl/LeyXjxkzzQTsuYibqO5k+7kS+E+S5pWJJF/cV/dR4CBJr1Am4Nws6fFlK0spXQtcDvxJft6OJTuPK7wfeK6kZ+e221MmTjuUTPvQJdMTLUo6nuxnhMqs4TPOA16UX19dsmv/iymlXZIeK+nxktpkNzc/YW9bG+dYA/V34B/R3vGIFwCklC4iExx8hUxE8NGBdi8k+53iG2S/e7xiiL7+hOxVyXfJBqj3rVSkLCTol4Gj8/pbyAaXlYHgVWRPaXeSDdwfqryn9/FWsovz02SDz3vIfm+8Nd+GV5I5ktcAv5xSuoXsPP8R2d3ibWSCjpWbiIuBrwE3SrpljO2ade4k0yp8UZkS/AtkQsRX5vXvIvv97MvA/wX+aaD968me9m4ns4UP9NW9l8w2bgCu5r5BtIzXkgmtvpC/Sv4MmTCRlNInyF6FXpwvc/EQ6ztr4Fq4Il/Xt8iEOp8Bvg0Mxvy/ikwUeRmZbfwla48Lax2r3yEbQK8mO17/k/y1JtmxezSZqO9jBW2HJqX0D8Cfk52LO8mepneklO4lG+iPJ7sW/wb4nfwpFbLrf1d+3F9KpgheeYo9D7gmfyU7yyr0ce25iP8O/JoylfVflyx3x4Ct/VFefiZwL5lzOhf4+5UGKaU7yURmzyX7ieTbwFOH2KYX5Pt5G9mN7Hv71nkd2Zuv15E56uvIbggbeX8vJxsrb8/Xc+EQ/UUU+oyUhTG/nuxpfzfZGPH8vM0WsvNwO9n4cCvw3/K69wBH5nZWGAU1DMp/UDfG9CHpCLKbsXZyzKupObbn+yd1fwI3xhhj9knswI0xxpga4lfoxhhjTA3xE7gxxhhTQ8Zy4JKeI+mbyvLADhOmYMzUsH2aWcb2acZl5FfoyiYn+BZZaMD1ZGEoJ6WUro7adNRNPRZWlasR3EdE5QDNam1Sozh+PywvnI12rXUFDYLyaPlUlmqg4rpoBOc3WF6K7aHRWC4sbwZ93PWtm25JKR1QWLnOTNI+R9yAaLsqLV+5HCCwz6hNZM/RelJJ3+F1U/GaqVpetc3i7bexdNfdoyaQGZuq9hmOne1gQslWPNHkcisYI1sVz1Fwrkc6R9GZqDyuxeNXZLbRuNYIxsJovGuqeD1lda2g/Mar7xhq7BxnOtHHAd9JKV0DIOmDZDF54QDZY4HH6+mryhtz8wVLgzZHGRRBC8VtljevmkgJgKWF4jwli/PFh2BxU+zB98wVW9We+WILWZorLl8s3gUWw7mXYGm+2HiWesXly3PFs1YqKO/O7Qn7nu8Vz1S4be4nheWXPP2tQ6cjXQcmZp8hjdhG1CyuiwZcddrFK2oX2626JXl3gnWlXnGb5bni5Zd7xeWLc/F+Ly4U1y1WvGb2LATXTMn91Z6F4mtgsaD8+295W8GSU6WSfUa22TrggQVLw/KBO8KO792veIz86bZi27x3c3DuNkXnLuyapeKuWZwLxrX5YgeXeoHTnYsj5Drd4rq57r2F5ZuC8i3d4vFuS7u4HGBb557i8taPC8v/8uh/GmrsHOcV+iHsnUv3eu7LNWzMRmP7NLOM7dOMzThP4EOhbO7TUwF6BI+cxmwQtk8zq9g2zVqM8wR+A3snyz80L9uLlNLZKaVjUkrHtFkrxbIxE8P2aWaZNe3TtmnWYhwHfhnwMEkPzifCeD7j5Zo1ZpLYPs0sY/s0YzPyK/SU0qKkl5ElzG8C56SUvlbWRo1GoWBNc8WqLbUDYQ+QAjFQqqpCb1ZU5lKisgyaVFZelmhjI41l1XVFavNSFfqGaXarM4p9hkSK8rIDUlEJjqKoiuoq9NB2KyraY2VwyX5H10DF8qrrKWtTWL7BtlzVPtVuFQrW0rbNhcsvBcJEgKVOsa0tByr05Ypq83IVevH4UnVMJRy/4r6rjnmRCj2ibPlGMHJX7WOQsX4DTyl9nGwaS2NmDtunmWVsn2ZcnInNGGOMqSF24MYYY0wNsQM3xhhjaogduDHGGFND1j2Ry140GoXpUSO1eZT6ESB1AxV6t1gymdqTyf+b1a1vXueyvkNVZpQDeELq9LK6KJ9vLSlIjxqqzYN0qVCSSjVqE/QR5k4vVcAHtl414mKkCI0JzRMwQoRG1SiQWtFqFaZHjdTmiwvx0L7UC1TogXA9BataDspLx84wV36kTq82l0OkTocytXnV5aPyeByM6srypw/D/cG0jTHGmH0OO3BjjDGmhtiBG2OMMTXEDtwYY4ypIXbgxhhjTA2xAzfGGGNqyHTDyJoNtLB6MpNwYpIgVCyrK455WG4XxymECfzbQQL/ILwsq4vKowlTipcPw8iacShEWBeG6QQJ/IPyRhS2ATQbxSEPZaFntUIqDvOKQrxKwsjCELOovBUYVVRe2nc0MUoUXlZxAqDSEMuovGLoZcVrpqxNGGJZI5ZbDe7db25VeTSuRaFiAItzxedisVdcvtQJxsiKYWdQcl7Dc1dcHI5fZROKBG2ica0VlWupeD2lE0EF2xtOTzUcfgI3xhhjaogduDHGGFND7MCNMcaYGmIHbowxxtQQO3BjjDGmhkx9MpPlzauVlClSxwYTk0CJ2rwbqc0DtWagQo/KIVZZhuUVlZehIrO0TaA2D8obgcKyTMUZqSwjtWbdkISKIiKCCUVGUqEHqnIF1wCtYHKeSGkOpKDv1Ko2yclyNJlJUF7eZjLLl14bVaI6ahY5kVrip9tW204U+RIpxKFEbR6Vd4vXE/Wx3I6P7XIrmLQkKCcc14LxKyiHeJyKVOjRRCOxOj3uux0q1z2ZiTHGGLPPYQdujDHG1BA7cGOMMaaG2IEbY4wxNWQsEZukXcCdwBKwmFI6ZhIbZcwksH2aWcb2acZlEir0p6aUbhlmwdQQSwudwvLC5QPlOJTlNo9yA08mzy/ECvUwR3qkqI0UmaX5nqNc6NXUms1gPa1msVoSoB3UdRqLYZsZYGj7REKdghOvKLd4rMaurDYv6hdIgQqdYP6ArO9AbV6xPLqWovkDoERtPrHIjRHmCShUoYermTZD2WdqwL2bV5+PeGyJ1xWNeaHafPWQnZcXH++ysTPOnx6cu1akNg/ympep0IO6duVc6NF64rEziu6J1OnD4lfoxhhjTA0Z14En4NOSrpB06iQ2yJgJYvs0s4zt04zFuK/Qj00p3SDpAcBFkr6RUrq0f4HcME8F6HW3jtmdMZWoZp+NhY3YRrPvUmqf/bbZWdi+UdtoZpixnsBTSjfk/28GLgAeV7DM2SmlY1JKx7RbHiDN9Khqnx2tzhJozHqxln3222ar57HTrGZkBy5pQdLmlc/As4CrJrVhxoyD7dPMMrZPMwnGeYV+IHCBsjzRLeADKaVPljVIDbE4v7rLKLdyCvL8QqyEDXObR8rLsDzsOlZlhrmBg/JItV6SSzhWa1bLhR6pNdslKs5OoLLslCjXN5DK9okE7YKTG6jNFeRIB0IVepjbvKLaPAVzAQAsd4rrovkDonza0fUXLV9WF9p6xciNaPnSuuia2Vgq2Wdqwp5Nq49tFLFSepyqjlOB2jwaI6PlAVI7GF8itXmwfLNVPOZESvOsrrhN1eiabrO4vFWiKO819hT3PaYKfWQHnlK6BnjUWL0bs07YPs0sY/s0k8BhZMYYY0wNsQM3xhhjaogduDHGGFND7MCNMcaYGjKJXOhDk5qwuGm1vDTMhV5yexGpXcM85YHCMlJSLnVLlLZVVeiBKjPMA12mmo3U5oGKM1abV1eU1zQX+vBIqFtwciO1eVku9GagNm8GRl1RbZ468aUbtQlzmwfzCoQRGqXzBBSXR1Ej4TUQ5syO+44Uzioqr9mjS2rAnoJQ8GiMLJ1PoWL0S+VxrUSFTtCHOtXGr1Yw3nUCdTpAN1KhV4yuiXKhd0vGwUht7lzoxhhjzD6IHbgxxhhTQ+zAjTHGmBpiB26MMcbUEDtwY4wxpobYgRtjjDE1ZLphZA2xZ271PUMcClE2aULQJigPw8vChPxh13GoTBA+EU7YMFIYRpTcPwgXawfhC0G4RbcVh0L0msUJ+eeC8trREHRWn9zQDhsl979BuFgKwstoRSFeQRhZyWQmS93iuqVuEC4WlEehX1E5VA+xjK6lpehaGiFEqfjamMkJTkJSA5YKZrtNjWAMic0jrFsOQlSjcxdOTFIyGZM6xeNOOH51isejTjBOReUQT0ISjWtReOxc897i9ZeEkXWDyUyi8mHxE7gxxhhTQ+zAjTHGmBpiB26MMcbUEDtwY4wxpobYgRtjjDE1ZMoqdNgzX6BgDUSto0xmEibqDxP4VyuHWAkbTYyy3I0UtYGKM1BkAjQC5WcrUJtHqsxeUB4pNQHmW8Xqy/uNCl0i9YafzCQ1SyYzCRTqKVCbR+XLFScmgRJVeTBBTzhpSVgedh22qRrtEU+gEV8bkcK5SMmsQL09qyTB4tzqbQ7HyJL9i1To0SRK4eRKwTgVTUwCo6jNg2iZdrVxDeIIm14w5kXjWqQ275UoynsK2qh4TB0WP4EbY4wxNcQO3BhjjKkhduDGGGNMDbEDN8YYY2rImg5c0jmSbpZ0VV/ZDkkXSfp2/n/7+m6mMcXYPs0sY/s068kwKvSdwFnAe/vKTgM+m1I6Q9Jp+ffXrrmmBizNrVapxrnQ41WFSsooz29UXjGvOZQo2gO1eZjbPMgZXKbijNTm7UCVGak151rFislIaQ6xKnO+MZ6Sckx2MiH7TA2xPFdgEJEKvSRXfwpzoUfRE4FyPIi2WO5UV6FXVpt3o/WX5EKP2lSO0AiumTKFcxCh0eusttuGpqZC38kk7LMBS/MF+xdG8JTsX2Q6zWpq8ygiptmMz9F6q82jcQ2qz+UQjndBLvSycTDKeV6mXB+GNZ/AU0qXArcNFJ8InJt/Phd43lhbYcyI2D7NLGP7NOvJqL+BH5hS2p1/vhE4cELbY8wksH2aWcb2aSbC2CK2lFKiZG4+SadKulzS5Yv33D1ud8ZUoop97ln88RS3zJhy++y3zaW77prylpk6MKoDv0nSQQD5/5ujBVNKZ6eUjkkpHdOaWxixO2MqMZJ9tlvzU9tAs08zlH3222Zz06apbqCpB6M68AuBk/PPJwP/PJnNMWYi2D7NLGP7NBNhTRW6pPOA44D9JV0PvAE4Azhf0ouBa4HfGKaz1IDFgoecNEIu9FC5XjHPb6Qoj3KqQ6yQDfM0R2rzbrHyshUoNQE6QV0vUGvOtwOFZaDWXChRoW9q/rRS+TSYpH3SEMu91Sr0yD4pyYUeKdSXgzZRbvMUqNAj5Xi2rvVVm0fLZ30E2xSozZeCayl1A+VzcM0AdHuBangDVegTs89GIvWKVOjBfpSMnVEeeAXq8UagTm8GyvFWyVwO0dwMk1Kbl0XRRGPbQqt4/IpU5VF5pDQHWGgU99HTeCr0NR14SumkoOrpY/VszASwfZpZxvZp1hNnYjPGGGNqiB24McYYU0PswI0xxpgaYgdujDHG1JBhcqFPjCRY7BVUjKJCD5SRoQo9WNdyoBCPVOtQkts8UF9Guc0jtXmnEyttJ6U239QuVkVubv0k7HtToNbc3Izb1IkksThXYECBorw0F3oYJRGp0IPySIUe5PAvW1ekEK+qNi9VofcCtXmU87xIWQ2oF+T878YRGvPdYnXwlu5q+2xOLxf6ZGgkGnOr9z1I0x+r04kV+I1AhR7lNm8F5VFe86xufdXmZVE0C0EO80nlPI+U5hCrzcvaDIOfwI0xxpgaYgdujDHG1BA7cGOMMaaG2IEbY4wxNcQO3BhjjKkhduDGGGNMDZlqGBkNWJpfHcIQT9Ycr6pyGFm0fBQuVhJGRjtI+h+Ut9rFYRVVJyYBWOgUhzBs6hSHI2xpF4d4ReVlE5Nsbd4TlN8/5nlPTVhcKDCgMMyxehhZNJlJZLcjhZFF4WLBZCbLUbhYtJ4gVKysbnkuaBOEkbWCcLH5XhwmtLVXbNM7uqvneW814gk3ZhEJOgXHREFIWFQO0AgmM4mOSRQu1moWj2vdoDyrCyYzicLLghCvcGKSIPQra1MtDHZTUD4fhH5F5bB+k5n4CdwYY4ypIXbgxhhjTA2xAzfGGGNqiB24McYYU0PswI0xxpgaMt3JTBrFKtUUiXkDtWRWF/QRqM3DdQVqcwUTkwA0g7pIbd4OVOVVJyaBWG0eTk4SqM2jSUu2toqV5gBbm6vVvABb7ieTmdAQi3OrDSuyz9LJdqIJUEK1ecXyYMKSrC4oD9XpwfLBBCTRxCQQq81TMDlJqxdMYjEXTEzSi5W+2wvU5gD7d+9a3a9ipfQs0mgsM1cwWUukNi8JkKAZqM2j8nZUHqjN242SyZgCVXkvUKdHE41EivJoeYjV5tHkJFUnLSmbmCRSqC9YhW6MMcbse9iBG2OMMTXEDtwYY4ypIXbgxhhjTA1Z04FLOkfSzZKu6is7XdINkq7M/05Y3800phjbp5llbJ9mPRlGhb4TOAt470D5mSmlN1fqrZFYnitQKEaKyRIlZaQqV6BCV5DPN1q+GSwP0I5ymwf5fLsV1eZzrViZGOUwj9Tm29pR/vKovFjJC7AtqNvWiNtMgZ1MyD5TA/bMrza6OEqifF1FhLnQK6rQo+UhzpMeqdBDtXknyGse5C8H4tzmkdp8vlidu3UuyGvei/PuP6BAbQ5wUOeHq8ra01Oh72QC9tlQYlOBCr0xQi70KOd5U8G5C5bvNIrPaackF3rUJlKPR+WhcrwkF3qc2zyIeGgUj5GxorwkD3tQNx8cj2FZ8wk8pXQpcNtYvRizTtg+zSxj+zTryTi/gb9M0lfyV0TbJ7ZFxkwG26eZZWyfZmxGdeDvAB4CHA3sBt4SLSjpVEmXS7p86c77x7STZuYZyT4X77F9mqkwlH3uZZs/jBMsmX2XkRx4SummlNJSSmkZeBfwuJJlz04pHZNSOqa5eWHU7TRmaEa1z9ac7dOsP8Pa5162uXVuuhtpasFIDlzSQX1ffwW4KlrWmGlj+zSzjO3TTIo1VeiSzgOOA/aXdD3wBuA4SUcDCdgFvGSo3hqgAhV6qJgsUaErUKE3IoVlRbV5lOcXoN0qrusFKvRIVR6VR3nNoUSFHuU2D9TmO1rFit2oHGBbs/gV87ZArTkNJmmfqQF7FgqMLsqFXmKfUc7zquXL0fKB0hxK1OOhOj3IX94NIjeCvOYArW613Oax2rw4suGBc3eGfT+wu1ptDnBoZ7WGrCxf9ySZlH02G4kt3eHnHIjU6RDngY/U5q1And4N8pdHywPMBSrxbqDGjlTlVfOXZ3XVcphHy29pRGr2klzoCubDKDlPw7CmA08pnVRQ/J6xejVmQtg+zSxj+zTriTOxGWOMMTXEDtwYY4ypIXbgxhhjTA2xAzfGGGNqyDC50CeGlOjOrVZeRyr0sny+jUiFHrRpBarydqBCL8vn241yngeqzPlWsTJyISiPFOUAm5pB7uhWtdzmkdp8v2aJCj1Qm29txLnb60RqwGJBKHioNi9ToQe3xrEKPVCOj5ALPVKVL3cCdXBQ3ugG10ygNAeY7wV5pXvFdhtb3hqZAAAFSElEQVTlNo/U5gd37wj7PrxzS2H5Ye1bV5V1AlXwrNLUcmEESjTeNUqU4M1ojAzaRIr9SM0eKcrL6nrBGBKpyrvB8pGiPFtXNRV6nL+8ePnNisfBhUDhP6+yCT/Wxk/gxhhjTA2xAzfGGGNqiB24McYYU0PswI0xxpgaYgdujDHG1BA7cGOMMaaGTDWMrNFYZr4gnKQRKOnLwsia0aQlQZtocpJOECJRNplJrxmEPARhYXPB8lFI2KZWHAoRTU4ShYttC8urT0yyIwjp2NEMYqNqRjaZSYH9jBJGFoWLheFlFcPI2nGYEO0gLLNTbNPNYF3dXmDn3XjCiK294hDI7d1iO3xAtzhsMZqYJAoVAzi8vXrSEoAjCkImu0xnMpNJ0dQy2zqrr80G1cPIotCzdhAWFi0fhX5F64E4/KsXhPVVDRfrlYZyVZu0JAwjC7Y1ChXL2hQPFvMqmZVoCPwEbowxxtQQO3BjjDGmhtiBG2OMMTXEDtwYY4ypIXbgxhhjTA2Zqgq92Uhsm4sn6hgkSq4PsUK9FSgBO0ES/WjSkmh5iFXlUXmUkD9SoW9uxsdoa6Ae3xK02dYIVOgjTEwSqc23NubCNnUim8xkeBV6pCgHIJhsJ1KbE5W3AkV5iQq92QompegEE0l0ArsNyrd0Y/vcEajN9w/U5gd1itXmh3aKFeVFE5OsUKQ2Bzi8tWlVWUe3h+uZRVpaZltr9bGNFOLNMhV6oFyP2kSq8qrlUKZCLx4jI6V7pDYvm8wkbBOUz0cTrwTHvGxikkhtPt/ohG2GwU/gxhhjTA2xAzfGGGNqiB24McYYU0PswI0xxpgaYgdujDHG1BClFOcbn3hn0g+Aa/Ov+wNxYuP1xX1PhwellA6YYn9jYfvcp/q2bY6G+54OQ9nnVB34Xh1Ll6eUjnHf+0bfdWNfPU/7at91Yl89R/tq32X4FboxxhhTQ+zAjTHGmBqykQ78bPe9T/VdN/bV87Sv9l0n9tVztK/2HbJhv4EbY4wxZnT8Ct0YY4ypIRviwCU9R9I3JX1H0mlT7nuXpK9KulLS5evc1zmSbpZ0VV/ZDkkXSfp2/n/7FPs+XdIN+b5fKemE9ei7ztg2bZuzjO3T9tnP1B24pCbwduB44EjgJElHTnkznppSOnoKYQE7gecMlJ0GfDal9DDgs/n3afUNcGa+70enlD6+Tn3XEtumbXOWsX3aPgfZiCfwxwHfSSldk1K6F/ggcOIGbMe6k1K6FBicF/FE4Nz887nA86bYtynHtmnbnGVsn7bPvdgIB34IcF3f9+vzsmmRgE9LukLSqVPsd4UDU0q78883AgdOuf+XSfpK/ppoXV5B1Rjbpm1zlrF92j73Yl8UsR2bUno02WuoP5T0lI3akJSFAEwzDOAdwEOAo4HdwFum2LdZG9umbXOWsX3OmH1uhAO/ATis7/uhedlUSCndkP+/GbiA7LXUNLlJ0kEA+f+bp9VxSummlNJSSmkZeBfT3/dZx7Zp25xlbJ+2z73YCAd+GfAwSQ+W1AGeD1w4jY4lLUjavPIZeBZwVXmriXMhcHL++WTgn6fV8Yrx5/wK09/3Wce2aducZWyfts+9aE27w5TSoqSXAZ8CmsA5KaWvTan7A4ELJEG27x9IKX1yvTqTdB5wHLC/pOuBNwBnAOdLejHZ7EK/McW+j5N0NNmrp13AS9aj77pi27RtzjK2T9vnIM7EZowxxtSQfVHEZowxxtQeO3BjjDGmhtiBG2OMMTXEDtwYY4ypIXbgxhhjTA2xAzfGGGNqiB24McYYU0PswI0xxpga8v8BLdNtfeLcgsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x216 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 20  # nb samples\n",
    "xs = np.zeros((n, 2))\n",
    "xs[:, 0] = np.arange(n) + 1\n",
    "xs[:, 1] = (np.arange(n) + 1) * -0.001  # to make it strictly convex...\n",
    "\n",
    "xt = np.zeros((n, 2))\n",
    "xt[:, 1] = np.arange(n) + 1\n",
    "\n",
    "a, b = ot.unif(n), ot.unif(n)  # uniform distribution on samples\n",
    "\n",
    "# loss matrix\n",
    "M1 = ot.dist(xs, xt, metric='euclidean')\n",
    "M1 /= M1.max()\n",
    "\n",
    "# loss matrix\n",
    "M2 = ot.dist(xs, xt, metric='sqeuclidean')\n",
    "M2 /= M2.max()\n",
    "\n",
    "# loss matrix\n",
    "Mp = np.sqrt(ot.dist(xs, xt, metric='euclidean'))\n",
    "Mp /= Mp.max()\n",
    "\n",
    "# Data\n",
    "pl.figure(1, figsize=(7, 3))\n",
    "pl.clf()\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "pl.title('Source and target distributions')\n",
    "\n",
    "\n",
    "# Cost matrices\n",
    "pl.figure(2, figsize=(7, 3))\n",
    "\n",
    "pl.subplot(1, 3, 1)\n",
    "pl.imshow(M1, interpolation='nearest')\n",
    "pl.title('Euclidean cost')\n",
    "\n",
    "pl.subplot(1, 3, 2)\n",
    "pl.imshow(M2, interpolation='nearest')\n",
    "pl.title('Squared Euclidean cost')\n",
    "\n",
    "pl.subplot(1, 3, 3)\n",
    "pl.imshow(Mp, interpolation='nearest')\n",
    "pl.title('Sqrt Euclidean cost')\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset 1 : Plot OT Matrices\n",
    "----------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAADQCAYAAADBEII/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXlYVNf5x79nBnBDcUNFRcF9STQquG/XpcYkxn3colDTxk6XpK0l1rYYA0lta/tLmjadkFXcvW7EbMYYx11Q3BU1LqhRARUVRWU/vz/ee5nLOAMzMMMwcL7PwwPM3HvmztzvnM89577nfRnnHEJCQkJCQkLeIZ2nD0BISEhISEjIcQlwCwkJCQkJeZEEuIWEhISEhLxIAtxCQkJCQkJeJAFuISEhISEhL5IAt5CQkJCQkBdJgNsDYowtY4y9pfw9hDF2zpFthYRcLcbYYsbYykp+TeF/IbeLMcYZYx2Uvz9gjEU7sq03qNqBmzEWyRg7yRh7xBhLZ4yZGGMNlec+YIxlKz95jLF8zf/f2GhrOGOsSLON+jPAVcfLOd/DOe/sqvaEHJcrvVIdJfzvPaqqXlY8dK2MbZYpx6X12HFXHgfn/Bec81hXtulJVStwM8bmA/g7gCgAAQD6A2gL4DvGmJ9y8vw55/4A/gpgnfo/53ysnWZvaLZRfw5UyhsScpvc5BWPijHm44Zmhf+ruKqql5304z+sPNbTXcdVHVRtwM0YawDgTQC/4Zxv5Zznc84vAzAACAHwkhte8zJjbJTm/xLTjoyxwYyx/Yyxe4yxHxljkTbaKHFFyhjrxRg7whh7wBhbB6C21fYvMMaOKW3uZ4z10Dz3R8bYRWXfFMbYRM1zkYyxvYyxfzLG7jLGUhljVRJA7pY7vMIYa8oY+1I5L3cYY3sYYzrluRLnlDG2VjNVHMkY22vVlnaK73nG2FHG2H3FQ4s124Uo277MGLsKYIfyeH+N744zxoZr9glljO1SjuU7AE2dfa+atoT/Pawq5uXhjLFrjLEFjLF0AGsAfAOgpWYk3dLJY3lixK71HWNMzxj7k+a8H2aMBdtop8QtF8ZYFGMsjTF2gzE212rbWopPrjLGMhjNWNRRnmukfDa3FB99yRhrrdl3J2MsljG2TzmebYyxcn/H7KnagBvAQNCXfJP2Qc55NoCvAYyuzINhjLUFmfY/AAIBPAPgWBn7+AFIALACQGMA6wFM1jzfC8CnAOYBaAIgDsAWxlgtZZOLAIaArrrfBLCSMRakeYl+AM6BOut/APiEMcYq9Ea9U+7wynwA10DnujmAPwHgZZ1TB/QQwBwADQE8D8DIGJtgtc0wAF0BjGGMtQLwFYC3lNf7A4CNjLFAZdvVAA6DPBALIMK5t+mYhP8rTVXNyy2U59qCfDsWJWdtbpTjeErT7wHMAPAcgAYA5gJ4VNoOjLFnQd+L0QA6AhhltcnfAHQCebYDgFYAFinP6QB8Bnp/bQA8BvBfq/1nAvgpgGYA/JTXcqmqE7ibArjNOS+w8Vwayj+yaKlceWp/6jmw30wA2znna5Sr4EzOeakdF2iKyxfAu8o+GwAc0jz/CoA4znkS57yQcx4PIFfZD5zz9ZzzG5zzIs75OgDnAfTV7H+Fc/4R57wQQDyAINAXs6bJHV7JB32ebZVzt4dTIYCyzmmp4pzv5JyfVM7pCdAoZpjVZos55w85549BI6yvOedfK/t8ByAZwHOMsTYAwgFEc85zOee7AXxRxiEI/1dtVTUvFwF4Q/HXYyde8w9WHot3cL+fAfgL5/wcJx3nnGeWsY8BwGec81Oc84cAFqtPKBdyrwD4Hef8Duf8Aej2wnQAUHy8kXP+SHnubTz5ffyMc/6D8v5l0AWAS1WdwH0bQFNm+75KkPJ8eXSDc97Q6uehA/sFg0YAzqglgOvKl0TVFc3fbQHM1xpceZ2WAMAYm6OZRrwH4CmU/OKmq39wztWrUn8nj7E6yB1eWQrgAoBtjLFLjLE/Ko+XdU5LFWOsH2PMrEzNZQH4BZ7sjH/U/N0WwFQrjwwGva+WAO5a+besYxH+r9qqal6+xTnPKcdr/tPKY47OBJXXZ9rvjPY9BAKoC+CwxkdblcfBGKvLGItjjF1hjN0HsBtAQ8aYXtNGuubvR3CDx6oTuA+Arr4naR9kjPmDpmu+d8NrPgSdZFUtNH//CKC9k+2lAWhlNX3XxqrNt60MXpdzvkaZmvwIwK8BNOGcNwRwCoC3TwW6Qy73Cuf8Aed8Pue8HYAXAfyeMTYSZZ/TEh5ijGk9BNDU9hYAwZzzAAAf4Mlzqu1IfwSwwsoj9Tjnf1OOpZHViLkNyi/hf8+rKnkZKOlFW/87K+vvhx4KRBWV12fa++Da93AbNP3dXeOxAE6BfQDdRugMoB/nvAGAoeqhOXkMFVK1ATfnPAt0X+s/jLFnGWO+jLEQ0FTFNdB9GVfrGIDpymuFAZiieW4VgFGMMQNjzIcx1oQxVtaUyQEABQBeVdqchJJTfR8B+IUyCmOMsXqMgpfqA6gH+pLcAgDG2E9BIw4hK7nDK4yCpjoonVoWgELQtGFZ5/Q4gO6MsWcYY7WhmbZTVB/AHc55DmOsL2gKujStBDCOMTaGUeBObUYBPq0551dA0+ZvMsb8GGODAYxz9r1qJPzvYVUxL9tSBoAmjLEAZ49D0Q8Aaivn2RfAXwDU0jz/MYBYxlhHxRM9GGNNymhTBhDJGOvGGKsL4A31Cc55Echn7zDGmgEAY6wVY2yMskl9ENjvMcYaa/etTFUbcAMA5/wfoECKfwK4DyAJdEU2knOeW85mtRGR6o8akBENutq7C/ryrNYcy1VQwMR8AHdAnVypSxw453mgK+dIZZ9p0ASdcM6TAfwcFAxxFzSdFak8lwLgX6AvVwaApwHsK+d7rvZyg1c6AtgOIBt0Dv7HOTc7cE5/ABCj7HsewN6SzeKXAGIYYw9AATJyGe/rRwDjlfd2S3lPUbB812eCgrTugDqd5WW8L+H/Kq6q4mU7x3YWFJdxSZl6thdV/rqVx24r+2eBvgMfA7gOGoFro8z/D/Sd2AZ6758AqFPGMX0D4F3QKowLym+tFiiPJyrT4dtBo2wo+9UBjcwTQdPolS5W8naFkJCQu8UYWwbgGuf8L54+FiGhikh42TOqViNuISEhISGh6i4BbiEhISEhIS+SmCoXEhISEhLyIokRt5CQkJCQkBfJHUUJ7Kpp06Y8JCSkMl9SyMM6fPjwbc55YNlbek7ClzVPwpdCVVGO+rJSwR0SEoLk5OTKfEkhD4sx5nCWME9J+LLmSfhSqCrKUV+KqXIhISEhISEvkgC3kJCQkJCQF0mAW0hISEhIyIskwC0kJCQkJORFEuAWEhISEhLyIglwCwkJCQkJeZEEuIWEhISEhLxIAtxCQkJCQkJeJAFuISEhISEhL5IAt5CQkJCQkBdJgFtISEhISMiLJMAtJCQkJCTkRRLgFhISEhIS8iIJcAsJCQkJCXmRBLiFhISEhIS8SALcQkJCQkJCXiQBbiEhISEhIS+SALeQkJCQkJAXSYBbSEhISEjIi1QmuBljwYwxM2MshTF2mjH2mvJ4Y8bYd4yx88rvRk698j/+AZjNJR8zm+lxIaEyJHwpVBUlfClUGXJkxF0AYD7nvBuA/gB+xRjrBuCPAL7nnHcE8L3yv+MKDwcMBosZzWb6PzzcqWaEaqyEL4WqooQvhdwvzrlTPwA+BzAawDkAQcpjQQDOlbVvnz59uFYF3+3guQFNee6CaM6bNuV8xw4uVL0EIJk76bHy/LjDl3dfFb6srvJGX/IdO3hew6Y86zXhy+oqR33p1D1uxlgIgF4AkgA055ynKU+lA2huZ59XGGPJjLHkW7dulXgus4eEg72N8Pt7LPJeNgKS5MzhCAkBcL0v2QgJB3oa0fC9WMAofClUPrnal7kDJRzpZ0SDf8fi/izhy5osh8HNGPMHsBHAbznn97XPKVcK3NZ+nPMPOedhnPOwwMDAEs81O23GgOMm7JWiUfi+CY+/NttqQkjIrtzhS90uM8KTTbg4MxowmZ68tygkVIbc4cta+8mXSWOi4fOxCbdk4cuaKofAzRjzBZlwFed8k/JwBmMsSHk+CMBNp15ZuUej3yAj+LMYbJomA9MMAt5CDsudvvxitoyU6TGALJe8tygkVIbc6Uvdehld1sVg609l1JtrEPCuoXIkqpwB+ATAGc75/2me2gIgQvk7AnQvx3EtXQosXAhIEtq2BQZHS9g3ZCEyopYiO9uploRqoNzty5vdJeTlgaYjFy6kx4WEylBl9JcBAcDItyQkj16I+9FLcf26Sw5dyIvkyIh7EIDZAEYwxo4pP88B+BuA0Yyx8wBGKf87rqgoYMkSFG6nK8a2l8yQkpZg34AoxMdDwFuoLLnVlyGpZgK32QwsWUKPCwmVrUrpLwOOmDF4zxIcGx2FFSsg4F3D5FPWBpzzvQCYnadHlvuVJQnZn8rQTzTgUYQRTdaZoN8gY3A7CatWAfHxQEQE4O9f7lcQqsZypy/zV8kYOdGAUylG4LCJpstFIJCQA3KnLwtWyyiYZMD1SUa0+coE3XoZo3pLuLYMWLECmD0baNWq3K8g5EXyaOY0/SgJKUONaPJ+LDKnGYunzWfNArKyIEbeQh6Rz2gJyeFG9P1WRJULVR3pRkpIHWNEm/hYXH3eWDxtHhkJ1KkDMfKuQfIouOskmtH7oAmHn4tG7c9MuLZCmTYX8BbyoNhOM8IO0WoHEVUuVFWk22VGl50mnJkSjSbrTTjxb2XaXMC7xslz4FaiJJkso9uGGOyYJ6PRLwwC3kKeleLLL+fIMEsiqlyoikjTX3ZeF4PDr8to/ycDTr5nG943bnj2cIXcK8+BWxMlWacOMOptCcefW4jHsUtx4QJt0rYtMHOmgLdQJUoTVV5UBBFVLlQ1pOkvdTpahXNh6kLU+s9S7N1LmwQEUFxQnTrA8uUC3tVZngO3EiV5e70ZnNO0+YCdS3D2hSisXYtieIeECHgLVaKiosCXLEHbSzSSKfpeRJULVQEp/eW9zbTaQbfLjB5fLUH6S1H4/nsUw7thQwHvmiDPgVuSkGmSUSfSgAszF4Er00Cj3pYQGAgBbyHPSJKQv1LGiDgDhu9YBEwziKhyIc9LkpC3UobfbANSplr6y8HREp56CgLeNUweDU5rPFlC+gQjOq6NxYVRRvDhNG0+Zw5Khffy5QLeQu6T3xgJR/sZMWx3LLJmiKhyoaohvzESHrxkxDNbYnF8oBF5g2jafOJElApvcc+7+smj4GY7zWi3jXJCt9xiQtLflGnzMuB9756At5AbZTajz0ETdg2Nhv9KEVUuVEVkNqP5RhMy5kWj43YTdkQr0+ZlwLt2bQHv6ibPgXvePGDiRDBZRruVMUh5Q8Yzb07Ej8/NE/AW8pwUX+75jYydI2Jw8i8y9Yrz5nn6yIRqshRfQpbR/IMYZLwnY/i/JyL1J/MEvGugPDriBtWmBWNAWBig9+G4eQvYuhXF8J49W8BbqJLFORo1oj+zslDsUyEhj0rjw3btAF9fjuxsYM0aCHjXMHkO3HFxKNqUgLyJBuT9cRHYNAN8tiQg8+04HDxogXfduo7B++FDj70ToeqkuDgUbEhAj7cpOK3/OwYgIQGIi/P0kQnVZMXFAQkJyJ9kQPbvF1FlxS0J8P00Dleu2If3vn20u4B39ZJHR9y3npKQ1MsIv7/HIv9nRrAREn7yE6B/fzgF77t3KdpcwFvIFWIjJBzsTcFpx/qL4DShqqHcgRIO9zXC/51YZM8mX/boAUyYALvw3r5dwLs6yqPgbp5ixsATJuwZHo2C/5qQ840ZjMFpeM+aJeAt5Drpd5sRlkzBaT33i+A0oaqhWvvN6HvYhMTR0dB9aMKdjeRLAe+aJ4+nPNVvkNHqkxhsNMjgBkMJePfrZxveTZsKeAu5SYovTy+i4LT1U0XKU6EqIMWXuvWU8vSrOTJqRxgcgnf37qXDOy3Ng+9LqFyqEilP27UDBv5Zwt4hC5ERtRSPHlHA2pgxtuE9Z46At5CbpPiywXiaHr8cKqEgSqQ8FfKwNP1lo0aUIvrgyIW4++elSE+nTezBe9Ik+/CuVYtihAS8vUseT3la9D1dMba7YsaIpCXY0y8Ky5fDKXhfvEhNinveQhWW4suQVPJlSKoZ7O8i5amQh2XVXzY6ZsbQfUtwZAT1l+WFd2SkgLc3yqMpTx9+JiN3ggF3X11UPG0+4E8SMjPhFLzXrLHAOzRUwFuoApIk5CyXUTSVosqnrKdpcxGgJuRRSRIK18jIm2jAtbmLiqfNR74lwdcXpcI7P1/Au7rJs+u4JQmnhxjR6D+xuKuklmzfHpg+HaXC+9tvHYe3WCom5KxqPSvhWH+KKk8OM+JcSwFtoSogScLFnxjR+rNYXBtH/WXjxjTlXRq8V68W8K5u8ii46x2k1JLJY6Ph94kJaatpGqgseCclOQ7vO3cEvIWcE9tpRq8kiioPSzbBd68ITBPyvPS7zei2y4TTk6PRaJ0JKe+TL8uC9+XL9uG9fz9tK+DtXSoT3IyxTxljNxljpzSPLWaMXWeMHVN+nnP6ldXC8OtldF0fg+0/l9Hg54YKw1t7z1vAu/rK3b7M/B9FlW+YKuMnH4uociHH5RZvqv2lLKPLuhgcnC+j7esGh+A9caJ9eH/3nYC3N8qREfcyAM/aePwdzvkzys/XTr+yJkqyXj2Kkjz67EJkv7kUqam0SXng3aSJgHcN0TK40ZfNp1PlpcuhEvYMXogiEVUu5LiWwdXe1PSXej0w9A0JP0xeCJ93lyIpiTYR8K45KhPcnPPdAO64/JWVKMm7m6giWL2DZgzatQSnx0Zh9WoIeAuVKnf7Mu9bM/z9Kap88N4lePCKiCoXckxu8abiy/ufm5GfT9Pmz3yzBNdnRGHrVjgFb+uAtW7dBLy9TRW5x/1rxtgJZVqokdN7SxLufECF4VNnK4Xh18sY/TatU3QU3n37Og7vGTMEvGuAKuzLx/EUVT7wW4oq3zBVxg+tRICaUIVVfm9KEvJWyvCZZcAZg9JfyjKGviGhSxeUCu+MDHpcveedmloS3pMn24a3WOdddVVecJsAtAfwDIA0AP+ytyFj7BXGWDJjLPnWrVslnms0ScKNF41otyoWqWOM4MNp2jwiAnbhfft2SXg/+6zj8G7XTsC7msslvqzznIQLo43o9y1FlV8OlXDpkluPW6j6yyFvluZLvzESsmYa0SMhFicHG1EwhKbNp0xBqfCOj7fAu2dPx+HdqJGAd1VVucDNOc/gnBdyzosAfASgbynbfsg5D+OchwUGBpZ4ju00o8N2E85Pj0bzzSYc/qcybW4F78uXafv27Qm8At5CtuQqX8JsRvfdlqjykFSzyOksVCE56s2yfBm02YS0n0ej/TYTzIvMKCiAgHcNVLnAzRgL0vw7EcApe9valVIYnskyOqyOwcm/yHgqeiKuvTDvCXivWuUcvLdtKwnvxo0FvGuCXO3LfWMoqnza2okYumqey45TqOapwt5UfAlZRtCHMbjxjowh/zcRl0bPcym8bd3z1sJb5DavGnJkOdgaAAcAdGaMXWOMvQzgH4yxk4yxEwAkAL8r16srheEZo0AznZ4jI4PqyDoK7xUrnoR3YmJJeEdEWOCtTnkKeHu33O1LvR4IDlb+BQcH8OCBSw5dqJrLbd5U+ksA6NgR8PHhuP+A+jVH4O3jU/Kety146/Ul4X3gAG2rwtvPT8C7KohxjRncrbCwMJ6cnFz8P99hRsEkA2A0wvdjE/g6GV8/lpCcDAwaBIwcSUB++JCuFu/epUIiISG0/4ULZNrAQKoaVrcueXvrVsqw1r8/VRljjOAeH0+gnjGDwA0QyNesIXPPmUMXC0KuE2PsMOc8zNPHUZqsfVm43YyiqQZcfc6IFgkmbJgq43KohOeeA8LDPXigQi6TN/oSZuovc+caUW+5CZBlHG0oYcsWSwyQjw9QWAhs2ACcPUuDmX79aPc7d4Bly+j5OXOA5s3p8ePHgYQESwCvry9ts2kTkJJCfeiAAbStmko6L4/63KAgCLlQjvrSo5nTMrpJSOxlhO/fYpH/MyPYCOocw8IoHZ+tkbf2nneHDmTWW7ecG3mvWSNG3kL2lTtQQlIvI9qvtgSnAcCZMx4+MKEardyBEpLDjaj3f7F4OIdSnvbqBbz4It0GtDfyPniQ9m/cmJZ56fXOjby3bbM/8laXmwlVrjwK7hZnzBhw3IQ9w6NR8F8TcrdSLW578J4zh5YpVATe6j3v0uD96JHHPhKhKqC6SWb0PVIyOA2ACFAT8qhq7SdfHhgVDRZnwt1N5Muy4P3NN87Be+1ax+BtvVZcqPLkOXArKfx8NspoHkdBQEVTDSXg3adPSXj7+zsG78ePLfAODy8Jb/UCoDR4x8cLeNdYKb4sWmMJTpuy3oCQVDNyc2mKUEio0qX4UrdeRsc1Mfhytoxacww24b1uXfnhPX489YmOwDsyUsDbU/IcuDUp/Dp1Avr9kVJLZkQtLQbv88+XD97Ll1vgPXasgLeQE1J8WXushB49KOVp0oiFGLCPUp6qKxOEhCpVmv6yaVNgRCz58u6flkJd7t2rFzBuHMX+OALviIgn4f3MMwLe3iDPgVtJ4cd30BVjp+tmjEhagp1hUVi5EgLeQp6Rxpfh4ZTytL95CQ4MopSnx497+PiEaqas+sumJ80Ytn8JkqUoxMejGN69e9uHd+fOJeHdpInz8O7aVcC7Kshz4FZSS+aMNyDrtUXF0+b9F0pIT4fD8A4IcB7e333nOLzFPe8aJsWXuRMMYG9QytP1U+TiALUrVzx8fEI1U5KEwjXkyxs/W1Q8bT4iVgJjcAjeU6dWHN6TJwt4VwV5NDitYIiEk4ONCHgvFlkzjcXT5tOmoRjeOTlPwnvHDgu8IyKch/eBA47B2zo/ulDNkN8Y8mWzOIoqv9LOkqc8JwfIzfXgwQnVWPHhEs6PMqLlJ7G4Md5YPG0eEQGXwvvmTXrcEXgnJtK2ImCtcuVRcNdPNiPskAkHn42G78cmZKxVps018F6x4kl4791bOfC2VdxEqPpLv9uM3kmWqPKnb5esxX34sIcOTKhGy2ePGU/tMeHkxGgErDHh3AfKtLmL4R0f7zi8v/3WAm97lcmEXC+PR5Xr1svoKsfg25dl1H/Z4BJ4q9OZAt5CTkvxpX6DjHu/p6jysZ9RVDljtMmJE549RKEaKMWXTKb+8sBvZbSeb8C5OAHvmqgqEVVevz4w6m0Jh8csxIM3luLqVdqkUyfAYHAe3qtWuQbepZUVFaqm0vhy9GjgSjsJ+4dTVLmaZPDmTaCoyLOHKVTDpPGljw8w/E0JZycsBPvn0uIZIBXegIB3dZfHo8rvf05XjPWTzRi0ewlOjonCqlUohnfnzq6Ft3ad99ixlOjFHrzLqgkuVA2l+DL7CzNq1QK6ppvRb8cSJA6OKh5xc265JSMkVCnS+LKggKbNe29bgquGKHz5JUrAOzKS/rYHb1kuP7xffFHAuyrIo1Hldz6QoZ9pwOU5lijJUW9L8PeHw/Du3dtxeE+bRibUwlvN0ibgLQQAkCQ8/EyGboYBVyIW4YUVBnw+k6LK69a1bLZvn+cOUagGSpKQt5J8edawCFyZNh/+poSOHeEUvM+ftw/vQ4doW3vwVhO9CHh7Vh4NTms0ScK1F4wIWRGLy2ONxdPmkZFwGN4vvOA4vDt2dB7ean50Ae+ao3ovkC/brYrFmWFG3H1GAudUeEHVlSslijUJCbldfmMk3J1uxFObY3F6iBGFQ2na3GCAU/B+4QXb8O7UCfj6awFvb5BHwc12mtFphwnnDNEI3GjCsXeUaXMreP/4I21fWfDevt12cRNreGtrggtVI5nN6Pi9CcnPRaOz2YTmKWY0bUp+U1VYaAlgFBKqFJnNaPW5CddfjkbotybsfMOMwkKUCm9b97z79LENb4OhfPBet64kvLt0KR3e6lpxofLLc+BWCsMzWUantTE4tlBGl4UTcX3cPAAl4b1yZeXCe/9+x+Ct1gQX8K5G0vgydHkMNhpkvPDJRLzwBflSp/nG7N3roWMUqnlSfAlZRquPY/DjP2UM+udEpP5kXgl4d+hQEt6Bga6Dt05nG95qfnQV3mqK1W+/fbImuK8vtSHgXTF5dMStzjUyBgwcCOj0HGlpwM6d9HT9+nSyS4O3NkmLFt5ms4C3UDml+LJJE6BlS4CDI0dJulJURJ0TIKbLhSpZGrN16QL4+HDcywLWr0cxvKdNKxvet2/T387COzLSeXhv3Srg7Q55DtxxceCbE5A/yYCCPy0Cm2aA7xcJuB4dh127LPBu0KB0eKel2Yb3nj0C3kLlUFwckEC+zF2wCEPfN0CenoC1w+MAkB/Ve92cizXdQpUkxZcFkw149AclRfQXCSj6XxzOnXMO3suWlQ1ve/e8Bbyrhjw64s7oJuFATyN8lsSi4OdGsBESXnyRlh04Cu+pUx2Ht5rb3J3wVmuCC3mv7veRkNTLiFr/iEXaeCMuh0po1Iiea9265LZ79lT+8QnVTOUOlHCojxF1/xWLRxEUzNu3Ly1rLQ3eR47Q/s7A28enYvC2rkxmC94+PuKed3nlUXC3OGPGwBMm7B4Wjfz/mJD3LWWncgbeXbo4Du/69V0Pb3U9rxqwpq4VF/D2XjU4bMaAY+TLZhtNCEk1Y9Qoeu7atZLbZmYCDx5U/jEK1TzV2m9Gv6Mm7B8ZDXxgQlYCBfOWBe8vvnAfvLW5zbXwXru2bHjbqgku5Jg8nvLUZ6OMpv+LwfopMgqnGFwO7169HIe3weA8vFetsp0fXcDbS6WmPN0oo9F/YiBPkTFlvQEtzlAnef8+bebra9lF9aeQkNukSRHdYXUMtsyS4fuSocLwLuuetzW8k5NpWzVgjTEBb0+oTHAzxj5ljN1kjJ3SPNaYMfYdY+y88ruR06+sSeHXrRsQFiVh96CFyIhaitxcMsS4cRWH97hx9uHdoEHJteJqilUBb++QW7yp8eXTT1Ougb2DFyLnraWoW5eW2ADkP1XHj4sgNSGL3O3LZs2AEbESEqWFyFy4FJmZtEl54M35k/B+/nn78P64xmhyAAAgAElEQVTqKwu8tfnRBbwrV46MuJcBeNbqsT8C+J5z3hHA98r/zsmqMHy3DDNGHlwCcx9KeZqbS1Mx7oR3RAS1tXKlY/Du04fgrdYE1yZpKQ3ejx87/ekIOaZlcLU3rXwZnm3G4L1L8F1PSnkaEAD4+QH37ll2KSy0dIhCQqgEXzY7bcbwA0twcGgU4uPhUniHhZUf3upyMwFv96pMcHPOdwO4Y/XweADxyt/xACY4/cqShJzlMnLHG3D/t4uKp83DoiRcuwa78N61i3avLHivXGmBt7YmuApvNejN+p63reImQq6VW7wpSShYLSN3ggGZv1qEZr8xYMNUGWyEhIcPyRPdutG59/Gx7Camy4VUucuXRWtl5E00IH2eJUW0FCOhsBBuhff69Y7DW7tWvFcv6ndLg7eaH10NWBPwdkzlvcfdnHOepvydDqC5vQ0ZY68wxpIZY8m31DOqKHeghOMDjWjw71jcn2UsnjafMgU24d2zJ3WQroD3zp2OwTsjw3F4N2wo4F0F5JA3S/Nl3iAJp4ca0eR/sbg+jqLKe/ak0faDBxZgFxRY9snOFoVHhEpVhX1ZOFTCOcmIFh/GIn0S9ZfNm1PfUxq8N2woCe/27Z2D9w8/PAnvjh0dg7eaH90evO0VNxHwLl0VDk7jnHMAdu/wcc4/5JyHcc7DAgMDSzwXcMSM8MMmJP0kGj4fm3BLVqbN7cD7xRddB+/du23D29Y9bwFv71Rp3izNl3WTzOidZMLx8dFosp6iygsKgB496Hm1w2vWrGSbX33l6ncgVB1VXl/67jXj6X0mnBgfjforTbjwMfWX1vC+o4z1VXifPVsS3tOnVxzeaorVisK7tMpkAt72VV5wZzDGggBA+X3T6RY0UZJd5Bh8Eymj3lyDS+CtLtkpD7zr1y8d3tZlRW3BOyBAwNuDqpg3FV8yWUb3DTE4OJ+iym/JZgQE0CbqbxXcarnP27ctATpCQlZymS+7ro/BvldlBL1msAnvZcucg/fRo7RtYCC14W54O1MTXBv0JmRRecG9BYCyoAARAD53ugVNlGRAADDqbQnJoxbifvRSXL9Om5QH3vXqUUCYO+Btrya4NbzVLG0C3h5Rxbyp8aWPDzD0DYoqb7dpKU6fpk1Gj6bzf/Ys+U2rTZsqevhC1VQu86WvLyDFSEgZvxBFf1+K48dpk/LCe8sWC7ybNXM/vC9ccA7e2kQvQiRHloOtAXAAQGfG2DXG2MsA/gZgNGPsPIBRyv/OSVMYHqBp88F7l+DY6CisWAHcuEGbOQvvyMjS4b1qVfnhPW2agHdVklu8qfEl54B+N/nywoSo4kx5eXlAcDB1gPn5Jdd0Z2SI0oU1Xe705aOvqCKY714zwr5bgtQpUUhIQJWBt3VlMgFv94jxSlyAGhYWxpPVyzIAdzaaUWuOAZlTjWjzlQmQZdzrJSE+nqA2Zw4VeQCAlBQyXOvWwKxZQK1aVPBhyxYy7fDhwLBhtO39+2Tahw+B2bMtaSrPniXDtWxJbdSuTebcsgU4dgwYOpTaYYyCkOLj6fesWUCbNtTGDz+Q4Vq0oLbVNlTTDhoEjBxJbWRnUxtZWcDMmUBICLVx4QLd61GnpurUqYQP30NijB3mnId5+jhKk7UvH31lBqYZcGWsEV12miBPllH3eQlNmwLbttEFWceOlo6usJAuFh8+pP8DA4Ff/tIDb0TIYXmjL/O+NaNgsgGXnzWi6y4TmCwjf7CENWso/fKECTSYAegCMj6ewBoZSVHbAC3B2rqVBjNTphAUCwqoP7p4kQZEvXrRtjdvUhs6HQFUzWFw6BAlY+nUiYDr40NtyDJFob/wAg1mAIJ+fDz1kRER9N0A6N76F1/QYGbaNGqjsJD653Pn6CKjb1/aNjOT2igspDasY0uqkxz1pUdTnjacKOHqc0a0iY/F1ecpSrJhQzo5derQqNTWyHv1atsj7927aVvtyNv6nveUKdSmOnrXZmmr6Mi7d28aedsqK2pv5C3WeVc91XlOwq3JRnTdEItDYUbc6CwhLw/o3588d/8+XQQWFtLFnrVu3RK1uoVcL78xEjKnGtFtYyzODDOicChNm8+YAYSG4omRd0QEAVUbsNavH/Dss46PvCMiaICkHXmHh1NeC3Xk7UhNcDHydq08Cm7dLjO67DThzBSK3j35Hk2blwbvyZMpctwa3j16UPyGNbzr1i0J765dLfBeudI+vAHn4V1WTfDVqy0pVtUrTW2iF6GqIbbTjLZfm3DvN9F4ao8JjY+bizucgACgbVvyDUAzP76+ltG2qs2bK/eYhWqAzGYEf2nCj5HRaPuNCXtilGlzK3irFetUeOfnuw7e6nIzLbxlufzwfuEFAe/yyHPgVgrDM1lG53UxOPy6jI6vT0Tai/MAlIS39p539+624T1+vG14R0S4D962AtbKgndZxU2EPCzFl5BlNHwvBkyWMW3dRIR9NA8JCXQui4qAl16izbdvp9smdevS/2qt7uxsS1EGIaEKS+PL4M9ikPo3Gf3/PhGXx8x7At6bN5cf3uo6b3vwXrasYvAGSsJbzY8u4O2cPDriVhM863TA4MGATs9x/QZBD7DAu3bt8sM7IKB88NamWLWGt3VNcAHvaiZN3EedOoCOcfj6AidP0r3De/co5sHfn3x3/TrlpFcD1tTlYd9+Sx2mkJBLpPHlU08BPj4cd+/RSoaiItfA29fXvfCOjKS/BbwrJs+BOy4OfHMC8icZUPjnRdBNN8BnSwKuLIzD9987D++8vMqBt62a4ALe1UhxcUBCAgomG5D3R0oteeD1BOyeFVccaXv/PgXXNG9O51Kt1a3+VvvXwkLq/ISEKiyNLx9HKSmiv0hA3ntxSEkBNm4sG95z5rgf3o7c8y4L3mp+dBXenToRvLVlRSMinqwJXpPk0RF3elcJB3oaof9rLApfMUI3UsLEiXQ16Sy8V62qHHhHRpYOb1s1wQW8vUsPwsiXfn+PxY0JRtzrRcFpISG0agCgiNisLOoA1ft3J05QQGStWpa2Tp0S2Z+EXKPcgRIO9jaizj9jkfNTCuYdOJDyCtiDt/aed4sWZcN740bH4W3rnrd1fnSDgeJ5SoN3aWVF9Xpqo7Sa4DUR3h4Fd9BZMwaeMGH3sGjk/8eE/G1m6HQoAe99+2hbb4G3rZrgAt7epfrJZgw6acKxF6MRsNqE3K1m5OTQc+qSmPBw6nAKC+mCLTiYvPnwIflHu7Z75UpR9lOo4qq134x+R03YNyIaRf8zFefAKA3eISGOw3vMGODMmSfh3a6dbXhb50e3B2+1uIm9e97Llgl4OyvPgVtJ4eezkYKA5MkyCqcYnoD39u0C3kKVKE0q3p4JMUh7V8YLyw1ofd6Mb7+1rLnv2JHWmgJ0X65WLTrfnTvTY/XrW5rMzga++65y34ZQNZPiS/0GGe1WxuDzmTJ0Mww24W19z7ttW8fg3b+/bXhPn26B97FjtG1p8C6tMtmXXz6ZHx0Q8HZWngO3JoVfjx7AM7+TsGvgQmRELS0G78SJBGVXwXvPHtrWGt5qitWuXakNV9zztgXvXr0I3mpZUW1ucwHvKiKNLxkDOvxcwo2IhRiwbykSEy3LvO7ftyS7aNCAOhmAzrGfH3WGOs2368CBmtOpCLlBGl8GBQHD35RwYPhC3Hp9aXFt+IEDgVGjgNOnS8J75kzXwfvzz8uGd3nKigKOTZur97ztwbum5Db3HLiVFH4w0xVjj0wzRh5cgh29o7BmjQW8kyY5B+9Jk+zDe8cO2/BesQIl8qOr8LaVpEULb3uVyezB215NcHvwVmuCC3hXoqx8CbMZbdcswYFBUZg+3bLsa/9+OicNGgCtWtGqCIBGJR07ku+Kiko2HR//5GNCQg7JypdBZ80YfmAJEodEYdkyFMN70CDvhbd1itU+fSw1wVV4a2uCW8NbjTWpCfD2HLglCbkrZOSMNyD794uKp82f+Z2EK1fgMLxr1SKwpSmVbp96yj68n37aOXhfv14xeJdVE7wseKtrxQW8K1GSBL6OfJn2yiJwgwFXl8q4HEopT195hUbU9+4B779PHdqtW5TmtlEj8srZs9Rhqqkj1frdjx6JxCxC5ZQkoWitjNwJBtz8xaLiafNhiyXk5qLawlvNj+4IvLWJXqo7vD0anPa4v4Rj/Y3wfycW2bONxdPmEybAYXhHRhK8ly8vG94TJjgPb21xExW8jsK7rMpktuCt1gQX8PacHvaVcHqoEUEfxeJgHyN+aCUBIA/o9ZT3uU0b6nwyMylq/PJl+l+ns+TGv3yZOrWCAkvbp04R2IWEnFXhUAlnhhnRLC4WNydTf9myJdVMcBTejtzzvnuXHvcmeNuqTLZ8uWW5WXWTZ3OVHzWj7xETEkdHQ/ehCXc2KtPmVvDOz68YvLXrvJ2Ft63KZO6Ed0SEBd7WWdpu3qQ2BLzdK/9DZvROMuH2L6PRY58JN9eRL1NSqANs0IDOwbRplkIIamBMfj4FCTVrRh2geq7UjGoAdT5qBysk5Kh895rR84AJx8ZFo94KE1I/JV86A28/v5LwPnmSttXCe9my8sNbrUxmD97WWdpcBW9HyopWJ3k8qly3XkantTH4ao6M2hEGm/Bevbpi8L561b3wtq4J7ip4r1z5JLwzMgS83SrFl0yW0fT9GNT+XMbMzw0ISTVj3z7gv/8lH92/T5s//TT97tLFcr4TEy3FR9TMaQ0bWl6Cc+Djj6kDExJySBpfdtsQg92/ktHsNwab8I6Pdxzemzc7B+9Nm0qHt7asqC1420qxKuDtvKpEVHnjxsDItyQcHLkQd/+8tLiesQrvy5ftw3v/ftrWlfBW86M7Am9bNcEFvL1YGl8CABshIff3FFXevz/56/Jl+uyPHKGgGIA6znnzaGR96hTBW6cj7zRoQB2YmlkNoPXeK1ZU/tsT8lJpfOnnR/3lqXELUfC3pTh9mjZR4Z2T4x54/+QnluVmroZ3+/YC3s7I41Hlj7+mK8bGx80Yum8JjoyIwvLlKAHviRPtw/u771wP78hIx+GtLStqD97asqKuhreaGETIRVJ8mbvVElVe998UVa4Gp4WH01NffEEj59q1aZq8WTPySa1a1FkVFVFnaTDQ9monqAarXblCtZGFhMqUpr8sKgL89pnR9/sluDAhChs3wil4b95cNrxnz34S3gMGuA/eamWyL76wJHopD7w7dqwZ8PZoVPndOBl8qgHXXl5UPG0+8i2qMWsP3tb3vLt1cx281RSrrob3ihW24b1qlW14W9cEt3fPW5sfXchFkiQUrKZkQMfHL0LBZANyl1NUuRqg2KkTbTpiBHV+OTk0jXj8OAWm5eYCs2ZRNrXcXGDtWoJ6vXrU0RUUWNZ4JyVZskkJCdmVJCFvpQwYDPhhOq12YDL1l8HBcBjeI0fSjFBZ8A4KssBbG7DmCLy1NcG18FYj1suCtzZLmyPwXr/eAm81P3p1h7dHg9MajJeQ+qwRrT+NxbUXjcXT5hERsAvv1NSS8J482XXw1uZHrwx4q2vFtfDWlhXVwtuRmuBCrlHRMAm3pxjRc0ss9j1txH9P07S52gGqWdEaN6YReOfO1PkkJFg8mJFhGWkXFdGI/OFDYMgQOtfa9dxffglcvFhJb07Ia+U3RsLNyUZ0WR+Lc5IRRcNo2ly9SHQE3oMHOw/vvDzn4J2QYBvey5aVH95qgR9reGuLm9QkeHsU3PrdZnTbbcLpSdFotNaEM/9Tps1rMLyta4ILeFe+/PaZ0eYrE/hfojH4pAlhD8iXhw6Rl1Qf3L9P50xNc/rcczRtDgDbtgGXLtEIvH794lvmMJtpKZm1tLdUhIRsymxG269NuDInGsFfmbDvLXMxeN0NbzXozRF4q8VNXAlvNWLdGt7WlclqCrwrBG7G2GXG2EnG2DHGWLJTOyuF4Zkso4scg4PzZYT+fiIyJswD8CS81QpLasCao/BWk7RUBN7OBKyVBu969Wzf87aXpU3Au3xyhS8hy2CxMdBvlCG9NxHjv56Hli2p4/riC9r01CnqhAID6f/69Wn3li2pA9m8mZ7PyCDP9epF51Jdo9+sWcmX/vRTi0eFqqfK7U2NL9vGx+DC2zLC/zoRV56dVwzemTPtw/vxY8fh3abNk/CeM8dxeGsrk1nDu6DAM/C2Lm7i7fB2xYhb4pw/wzkPc3pPpWSSXk/LZ3Q6jh+v0X0/oCS84+Mt8O7Z0za8bd3zbtSo4vBWLwBUeKspVm0laSkN3pGRTxY30cLbVn50Ae9yq8K+1P6v9yE/vvoqjW58femc/ec/tLoBIH8yRr4oLKRsfWpRkk8+sVQWe/pp6iTVgDbty370kSgDWgNUPm9qfNmzJ6D34ci8Q3ArKqI+zh6858xxHN4zZ7oP3hERtuFtXRO8LHjbKytqD962KpN5M7w9N1WuKQxf9JdF0M8wQL8lARej4rB1a/ngrdc7D++JE23D215N8Dp1qA3r4ibXrlEb5YV3aZXJHIW3CFhzgTS+zHmdgiaRkID9c+KQl0fnpUMH6gxbtqQpcLUj3LuXIK6u2a5dG/jVr+i85+eTJ3186H72z39OPrl5k34zRvtwDnz4oRh5C1lJ8WXhZANyF5Avfb9IwMN/xeH4cffC+9Qp2tZReNuqCW4L3vZqgjsCb7Um+LJlFYM34J3wrii4OYBtjLHDjLFXbG3AGHuFMZbMGEu+ZfXppHWRsO9pI3Rvx6JonhH6URKmTCGQ2YK3j4/z8D5wgLbVwlub2/zpp23DW1sT3FF4q1nabMFbLStqD95llRV1BN5qcRMB74r58lE/Cft7GFF7aSwO9zNit55uUGs/z/r1qRMcOpRG4UFB5L39+6njAaijKSqi2zv5+TQCr1uXcpavWEEXAAB5k3NLDe+iIhp5qzkAhKqVSvVmab7MHSghqbcRtf4Ri9y5FMw7bBgwfDhswrt1a9fAe9Mmx+Btrya4PXjbqgluC97qOm93wDsykv72NnhXFNyDOee9AYwF8CvG2FDrDTjnH3LOwzjnYYHqzUBFQWfNGHjChF1Do5H3ngmF283Q62EX3pGRdBK097zLgve2bU/C28/PeXhbFzdxBt7WNcErCm9AwLsMVciXdZPMGHzKhB8jo9FtF6WWTE+n+IZt2+iz9ven4DTO6TyHhFBH8+qr1JHq9dR5vPMO8OABtcsY8Otfk/9yc+l5gEbmISGWLGt0fHTPW+Q1r3Yq1Zul+bLWfjP6HTVhrxSNovdNePglBU3ag/esWRZ4p6RQG9b3vLOylIPSwDshoXzw1tYEdyW8fX0FvK1VIXBzzq8rv28C2Aygr8M7Kyn8fDfJaPBuDNZNkpE/yWAT3gcP0i4qvPV6+/Beu7bi8NYWN1HhbasymfU9b3vw1tYE18JbWxO8IvD29yd4W5cVTU+vmUlaXOFL3XoZwZ/FoM4WGbO/NKD3fTN0OrqQ/OwzyvBUWEijmYICClArKKBOZNgwyjTFGHWUajGH776jdae9epG/xo2j83/9OnV+tWrRj1br1llu+Qh5v8rtTcWX+g0yQpbHYNN0GbrpBofhvWGDBd6tWlngvWxZSXiPGEH3tZ2Ft7YmuApve9Pm6vdBDVgrDd7a/Ohlwbu0e97q6N0ReKtBb1VZ5QY3Y6weY6y++jeAnwA45XADmhR+vXoBPV6TsGvAQqTPX1pcNF2F9zffOA7vS5cqDm/rymSlwdtWTXAtvK1rgmvhrc2P7gy8tcVN6tenz8ORmuA1Qa70JQBAkqD700L03bUUtWpRAqvJky1BZRs30i5qEI96Dtu2pVHzgAHAb39LAH/4kDqRkyfpudu3AaORfJqTQ+c5N9dSXUzVd9+JcqDVQRXypsaXrVsDwxZL2Dd0IW69vrQ4b74r4D1kSPngrS1uosLbVmWy0FDyclllRVV4Wxc3mTaN1orbgrd1lrbwcIpY/+EHyrBWFrzVe97LllV9eFdkxN0cwF7G2HEABwF8xTl3PIGjVWH4XvfMGHFwCbb3isLatSgT3hERroW3ulbc1fC2rgku4O12udSXMJuBJUtw1RCFvDw61089RcE4AHWUTz1lqf2bkECftfrFv3qVZldGjKD/hw61gPnAAdq+Y0c6v2PGkB+vXaOlg1qdOAF88AF5SchrVX5vWvmy9XkzpKQl2DcwCvHxqFHwVhO9bNnieFlRR+CtZmkDqj68yw1uzvklznlP5ac75/xtpxpQUvjljDfg4fxFxdPmT78q4eJFmiIsDd5NmtiH9/jxtuHdtatteFsnevE0vNW14gLezssVvuTrZOROMOBq5CIUTjGgcI2M7HAJeXmWFTlq9rT69WnKe/58On9NmtC097Zt9PyBAxTc6O9PHej9+9R5jR1Lz6elWe5jHzlCz+l0JWt4qyVBMzKAf/3LcpEg5F2qkDc1vrz9y0XF0+ZD35CQnQ2vg3dIiGvg7UxNcFvw7tDBO+GtX7x4caW92Icffrj4lVcsgZRZjUNx6kA22q2KxcNfzIffvLkICqIRSmIidWrdutGH3LUrdViJidRBtmpFvzt1opN/7BidBH9/MkHDhrTt9euWNrp0oeCDpCQycnAw3afu0oUMevQo3UPx96f7L40bUxs//liyjTt36HFfXzJx7dr0+OnT1Pm2a0cderNm1JEnJdFFQPfu1EbnzmTupCTqpNu2tbSRkkImCg2lzyEwkNpJSqJc7epxdO5MXzI1eC8khN5T165k8uRkeiwggKaBWrSgbVNTLW1Uht588820xYsXf1g5r1Y+Wfsyp0UoTh7IRrcNsdgTPh+ras3F3bvUMQUG0rmtW5cuvIKC6HNmjM5PQQFNf3fvThdc9+7RRWRyMn3mGRl0bjt2pPPcpg0Vf7h+nTqa48fpgjI7m+6Tp6VRR6PT0UVDYSG1pX4HhMonb/RlfqtQnE7KRvvVsbgdMR91fz0XAQHkv+Rk+t537Ur9gOrJpCQCeqdO1F9160Z9UWIieTkwkPqZ0FDyY0oK9UO1a1O/pNdTG3fvUp/j40NtXL1KbTRpQv1T/frUdx45QkBX2wgOJuAnJpK/u3ShNrp3p341KYn83rw59bvt21M/fPIkvV6dOnSxUbs2tXHrlqWNbt3oe5OYSP19ixY0U9WxI7Vx4gS9by0vkpLoO9i1q6WNtDRqo0ED+j6rbRw/Tj9qG5UhR33pUXDXSTSj1XsLkDhoPlokmPCwazhqdw0tFd4ZGZUL70aNaFtreGdm0uOVBe/AwCfh3akTfSkdhXfz5pUPb2/sIH33mtHq3wuQ/+p8tPnahHrDwnHNJxQ5OXRu9u2jUXJuLsE0NJTOXXo6cOECBfn4+1OHlZJiWZpz5w5FmB89Sl7196cOcNQoyrt88CCdp9xcmqG5fp1mZ9S13tr85hcu0Hns3t0yIhdyXN7oS/1uM1q8swBHR1B/md46HAHPhJYKb6BqwNvXl9qoCvBOTHQc3seOVS68qz641cLw62XU+eVcfJsZjq6LDTUW3nq98/D29a368Pa6DlLxJWQZ+p/Nha5vOFr9zoDAseE48SAUEyaQr3JyCMT37pG/jh2j6casLOpkGjcmTyUmkkf79SMI799PvmjUiO5lFxVRDvSsLDrPV68CL79MvklNpYsBzmm7Fi0owE1VVha1HxRkqQsu5Ji81ZdMltHo93NhfhCObm8aqiS827V7Et5t2gh4O6KqD+5f/YpuJsycibp1geb9Q3HkhA9qyStRMG0W6tVDheDdsSN98I7Au3ZtMoWj8FbBq4W3nx9dAKjgPXXKOXgnJlZPeHtdB6nxJQD68H184CevxIHQWRgyhGIgevSwxCIMHUpwTUujqfKTJ2n0nJFBoM3JoX18fWmf27eBuXMp4vzcORpdZ2VZYizOnSP/6vXkr7Aw8ur9+3ReObfcay8qotdLTycvidG3Y/JmX/r4AG2GhuLEaR/4rV+JW6NnoWlTlBveSUkVg7faRtOm7oN3ly41A95VH9wtWwKvvYbcHuHw6RCKuklmtP2/17D9+Xex+2poMXjLC++jRx2Hd2JiSXh37uw8vBMT7cO7ffsn4X31asn71XfvOjfyTkykL4wW3s7c864seHtdB6n4sqhPOFi7UBrpvPYaMv/yLg7fCcXTT1tSml65QsCcNo380KcPTaN36UKeS0+nqfF792ikfekSnd8bN8gTjRqR306donz1ffqQpzMzKZDm9m26T/n4MV1HXLhA0+jqlLlebwF4ZiYFwjVrZsmJLmRf3urL3B7h0LcPpds5S19D4rR38f3FUDRvDqfgrQXv5cuOw1vNZaCFd/fultG7LXifPm0b3nfulARveeB9+3b1gnfVB3doKO51DIduhgG3UrNRP3YB2HoZzadLxR9QdYZ3YmLNgLfXdZChoSjqE46cFw04tjcbTf6xACejZaR3lZCaSp+xmtAqLY0SqgweTB2ajw/5pXFjWoHQty+dl0uXaIR+/74lw93RozSyBqjzyc6mSN4uXWi03q4djdLv3KFzeuQI+Tgnh85d27ZPRpcXFVEnee4cHad1Mhchi7zRl3k9w1E01YDUE9lo/I8FYLKMNhHky6QkOAXvrKyS8E5NJd9VJrwTEysP3h06EHirOryrPrgB+HQMxQ9HstFhdSxuzJyP+q/OLQavgHf1gLfXdZAA8luHIu18NrpvisURaT6+bj4Xqan03Jkz1Aldvkyfc2YmdUQBAdShXbpEI+wwpe5T7dr0mffrBzz/PEWKnztHI+mGDWlknZ9PHdjx49Re3bpUiOT552lJztmztE2jRvSajx7RSD44mF6rbt2S6VKzs8kXWVnkfZ3nSglVWXmjL3XtQ3HtTDY6rInF+XHz0fgPc4unvFV4t2hRPnh3726Bd7Nm1Q/e/v7eAW+vALdulxnN/rUAp8bMR1CCCalNw9E0LLRC8K5Xr3Lhfe3akwFr7oL36dPOw5sxz8LbGztI/W4zGi1ZAMyfj1ZbTBj8u6JnhcYAABn3SURBVHAEDw3FyZP05a9Xj0a7arKckydprfbp0zStnZlJn6tOR6PvgwfJD2oHl5dH8I6IIDA3b04do78/+VgdSR85QtPtQUF0fjt0AF56iTo2NUKdc4K2On2vXf+dnk5T9zodeVGtQCbknb5kO81ouGQBLk+aj6DPTUhGOFoPCXUY3t26OQbvpCTn4H3vnvvg3bhxxeHdqJFteKtBb1UJ3lUf3PPmAdHRYBs3oskf5mJvTjh6LJqIu0dSUW/auFLhXb9+xeEdEEDb3rhhOVFdulCnaR2wVhq8ExMrD95du7oW3mobWnhrl5u5Ql7XQSq+xMaNFEEWHg7d5IkIuJOK3QHj0KcPJXPo14/8dOQIpVYMCaGp6jt3CMwpKRQtvm8fff63bhFkHz2iz/zkSfJPcDCdvytX6Lnf/tYyRZ6ZSW1dvEiHlpZG5yc01JIoY+xY2vbOHYK2CmcfHzoezqkzPnCAPNe8uQA44L2+ZBs3ouHv5uJU7XD0eGMiru9NRcBL4xyCd0pK+eEdEuIZeCcmloR3u3bUD586VTLavFYt2/C+ds15eN+8WRLeN27Q4wEB7oe3o7707CSaElnj40OpoXU6jitXy67c0qcP8MILFKwjy5YMa1On0on4+mvqNAFLhjWdjtpQRzO9elEGMjVLm5phTc3S9u23dLKAkjXBtRnWevR4sjKZTme7JnjDhvRebNUEnzTJdk3wp58GduwA9uyhbQMCqA21uIka1dytGx33tWuUpa2smuBqYZIVK0rWBJ86lUy6alX1z7BWqtSIL83/TPmmaFOONm5Mvxs1orXYM2bQD0AZoyZMIMA3aED7mc3k11WraJvvv6frg927qVO5d4/82LQp3SP39aXO8Xe/o4x7vr7k3717qbN9/Jh82r49XQAAtL2vb8mRt3rcn38O/P3v9P2yfotCXiDlpDFG/Zfeh+N2JmX+4pxA+NJL5CVZtsRQBAfT49nZlA1MrVY3fDhlWTt2jDKscU7900svUSzchg10kQ8QHF96iS4utVXFhg6lvvvECfKXmmFt1izbNcG1lcnUDGuDBtkuKzpjBvl582ZLLYCgIGojL6/0muBqhjV7lcnmzKFtli2zZGnr29fxmuCBgdQG554pTOK5Efe4cUD//iiYbAAeZEO/cAGwcSN2Pf0bJCbSFVrLlnQlo81io46aW7a0jLzT0y1XSPZG3p070/7apWJBQbZH3mqWNkdH3o0b04hGO/Lu2pVOZmlLxZwZefv4WNaKe9PI2+tGNhpfPkzPhu9fFoBv2Ajda7/B/v10rtQ62j4+NKJu2pTOJUCf7d699Dn37UuPBwWRZyZPpkC24GDqRO/dowukc+csF1CnTlHHkZZG5+/8eXq9rl3pdU6doun1AQNo31u3yLtqR3rvHnV29etTMJyvb8nELYWFdF993z6Ce5s2NfMeuLf6snCKAQV3s+HzpwXQbdqIs6N+g6QkgrG6zKt7d4q1qMoj78OHnxx5+/jYjzbXjrzVNqxH3sHBFR95q329oyPvDh1cm2Gt6k+VA0irHYpje7IRujIWRb+bD/3P5pa4r6DC214KutLgnZ7uOnhrT6gW3mqiF0fgrSZ6KS+8ExPtw1tdK14eeKeklA5vNfiuvPK6DhJAQXAo9m3NRqd1sdgdPh8rfefiyBGayXjwgM7HzZv0mV65QjM1XbrQbzWy3NeXPl+AfLZ/P32+zzxDHUfz5nT+x46lkXmXLtTenTt0DtLTLbNDJ04QaDMz6fxduEDnMDycRi8ZGVSgJDiYvju3b1vyVhcVkV/q1y85i1JURD7bs4deR+30aoq80Ze5LUNxaEc2QlbEIu838+Hz87kIDaWLMWfh7UzAWmnwVtvwFLzVteLVBd5eAe76yZTydG9fCrZg4bSmu6z8sceOlQ3vbt3KD29b+dFtwfvIEcfhnZhYEt6dO9uGd+PG9td5W8NbG7BmC95qGyq81SxttuB9+DA9Zg1vNa1meeHtjR0kzGa0eX8B7kTOR8fvTQgYGQ5du1Dcvk2dZHo6jYTPnKEpuzt3CIAHD1IHoGZVU597/Jh8dO8e3eZhjM7VkSO0Tc+e9Lm3akVt9O5Na8MHDKCO9OZNum3COY3Uc3NplK52cDodnadWraittDQC86BBlvvu1rc+tGvA1aDMkyfpuJo0qf73wb3Rlz57zGj53gLs708pTwueCYdvp9AS8M7Opn6tNHi3bWsb3pw7D+/Tp0vCW82ProW3NsOa2j85Am9bWdq0AWvWKVbtwVvLBHvwtpWlzVl4ay8Aygvvqg9uTWH4m8/Pxc5sSnlaGfC2lWHNVfC2F7BmDW+1DWt4qxcAFYV306ZPwlubYtUa3mobWnirxU0qAm+v6yDNZrBpBujWy6j767nQ9wtH0GsGdJkdjjOPKdGF0UhQ7dGDoJeXR2uw1SVhDx7QvcArV2ha+uRJ6lCzs2ka/dgx6jQLC8knPj60j05H7V2+TPfG1fN86BDBdNo0qrikjpb79SPfck7n9epVSw71oiJqu2VL8tvNm3SuQ0NpNF5Y+OTn8PixJRd7ejqdf+vyotVF3uhLGMiX/KeWFNGVAW9tkpay4B0S8iS81cA5W/BWI9ZtwTsxsXLhfeRIxeCtDVgrL7yrPrg1KfxatQLyWobi7Hkf1Nu8EnV/PqvMcPyqCu8WLaoGvNUMWqXBWy1uYg/e2spk5YW313WQdlKeYuVKnOwxC4WFNN3t42NZFnb1KgXudOpEwYZ16tCI+Je/pFFvt260rXoe1GC17Gz6fekSAfPoUYJqbi595ikp1GGq68P1erov3bIl/X/7NsE8LIyCH8+do2MLC6NzevMmbZOWRh1yQQGNvps0oe/Ao0e0nZ9fyWA2zmm/Q4csnXjz5uS16iJv9mXDhkBAz1AcO+mDWutXwidiVvFFWUGB6+Gtgtcd8NYuN1PbqMnwrvrgVlL4ITwcCA1Fqx8o5ennw9/FubzQcpVdcwTe/v6WoDcV3sePW9bnlgfetoqbqPC2XivuCnjfufNkwJq94iaehrfXdZBWvlRTnuLdd3HmcSgeP6bpblW3btG0uZolDaCRz9GjdB5atyY/tWxJI9nu3em+9jPP0Ig5MZE+90mTyI/BwfQ5161LML5/n34KCujxEyfo/OTk0E9yMkWiP3pEHe7Fi/RaPXqQH8+cIT9OmECeTEujTr2oiNosLLSsylCXkKklRAF6Tu2gEhPpPnvjxnR83jyd7u2+bHjUjLbvvIZvRr+LpIzQ4mAzR+CtFqWpCLzVLG0qvNWgN1vwzspyvLiJJ+CtXW7maXhXfXCHhiL/mXAUTDYgL5Oid/UbZOQNkiq0EL5lSzoBpUWbW8NbvQBwB7wPHHAO3mqHbw1vWzXBywPv0sqK1qpVMkubK+DtdR1kaCh4WDjyJhpw+WQ2/N9agDNvyrjSTsK1awTRoCCaVs7Npf/Pn6d70P7+1IQaWd68OX2uAH1eZ88SbHv0oMfUqfFLl2g5WdOmNCOUk0OPRUbS1PigQfTY9eu0TKxbN/KeCvXatclT9+5Rh3v5smW9LEDT8OfOUccZFET/5+RQ29270zE8fkzeZ+zJpWSq1Pv7ycl0T//SJXosMND7ipt4qy/zJxmQdS0bdRYvgG69jAbjJSQnk7fcBW/rwiTWKVZV8JZWE9yd8Ha0Jrh1cRMV3gcOlB/eXbuWDm9tljZHVPXBDeBuQChOHchGu1WxeGScD995c12SxaY6wjsx0TXwdrYmeEXh7XUdJIDC4FAc2ZWNHp/HYm+/+fimxVycP0/Ay8sjnxw5QlPJ58/TPocP03k+eJCez80lr125QttcvEj737hBn9utW3T+8vOpM23YkD77oiLyTHIyebNNG2o/KIheT6ej9bfBweSFo0dp+1deoTW1ISHkocBA+r9FC2rz7l16/du3CdIA+UktXKLXE8z1empDp6Pt9HrypjalKkCdelYWvbc9eywXl76+5P2KLiN0t7zRl/mtQnFyfzbar47FnZ/OR51fzUXDhuSF6g5vV9QE9/OzzBqVleiltLKiSUmOw9u6uElZ8gpw102iqPIDAylK8lG3cNTuEupxeFfknrcr4W29Vrw0eDtSE9wT8PbGDlK3y4zW/6GUp22/MWHQa+HoOy0Ujx7Rl372bPJAp07kEbX4SLt25Jc6daizVJP63L1LnlODwi5doqC1lBQ6n4AlSnz/fksCosuXyQdHjlCnzDl1VBkZ9Ds9nT7/ixct0996PZ3vM2eoMxsyhO556/UE6fBwSo7RvbulelloKHmNMcuSNBXuakpVgJ5Xp8et134XFNBnc/o03RI4cIDeY1YWfR7qaL6qyBt9qd9tRtC7C3BkOPWXt9qGo0GP0BoBb2drgjsDb7WsaFWAd9UHt1oYfr0Mv19QlGS3Nw143D0ctazgbX1fwd3wdjRgLSODTqCarN6V8LaX6MUevBMTKx/ely+XnaTF6zpIxZeQZWDuXLDwcOhnGOA3KBxptUORmkpZzNQlc61bU+azp56i6e7OnekzefiQzvVrr1FhkYED6b52YiJliRo/nrJf9exJPioqAp57js6VmgYyM5N81aABHRpj1AFnZdGIPTWVtgHoO3H6NHlPrUB2/TpBdPdumhLU6ei548fpOXWEn55OFxQdOtBr379Po+/QULp3r9NZpuEDA2m/nJzSs68VFlI7V66Qf3btoveekkJ+ZYw6UT8/953X0uStvmSyjEa/n4vvsyiq3Ba8z52zQLMseB88SP2TFt6HDtEFQFWAt6M1wV0F79KKm9iCt7YmeNeujpcVtadKSXnKGHuWMXaOMXaBMfZHp3Y+dIg6R0lC8+aAFCNhyywZxz46hDt3gMWLqdMYO7ZkCrq33qJI2vbtgS++sKSge/99CrrUpqBbvJhGG88/T1N6anrUt96i9J6dOgFffUUmA4D//pfaYMySYnXxYupgx42zpEctKABiYy0pVr/5hr4AAPDee9SGXk+pTTMyqI2ePamzvnQJWLuWRjGxsRSU1LUrsG0bnXAA+Pe/6f6mNsXq4sV0H3XCBDL6mjU0bRsTQ/c9n3oK2L6dOmoAePddOo7atSm16Y0b1Eb37pTB68cfKfWm2sb48XTv1Wymjh4A3nmH2qhbF1i50pJidd06SrF6/bolxWolXv+VKVf5EgD9lmXg0CH4+dHno017qtfTlzMurmQzgYHkk4ULLY81aECd47//Tb+bNaNOoHdvAv3KleT5IUPIW3XqUKeXkkKj/FdeoQsAdXQbHQ386U/0HEAd2N27FBA/fjx1PH5+tM/Zs+TBgADqwB8+pJ+tW8lnGRk05b1nD6WYBKhT3raNvjtFRfTe1Sl+vZ68pdPR4z4+5BM1QM9sfnKEnZtLaX337ydP/vOfwJtvUhrWF14A1q8nwF+4QIB44w3bp8ie12w9XpV8CVTAmxpf1qkDjHpbwo55Ms6uOIQLF2iTZcvo3N+7R/1GdjZ9viNHUpzE4cOUDppz4G9/o5UQzZtTsz/8QG18+ik9fv8+9YEPHlAbw4eTL48epX5X20ZQEJ27s2epjY8/tqRYjY+nthYvpvSqw4bRheOWLeSpJUtoBqh1a0qPmpJCbXz0Eflam2J18WI6hhEjCKQJCdTGX/9K77tNG0p3euoUtREXR2lJc3Pps7l3j9oYOJBSEp8+TdsXFQFvv01ttG1L7Z44QW188AG1kZ9Px6GyqX9/Snp05gwdd2EhtTF9Ol18f/45Dd4A4H//o360sJDayMysuC8ZL2fSYsaYHsAPAEYDuAbgEIAZnPMUe/uEhYXxZJWSNpSRQYbT64H58y1X9AcPEhy7dKEPRl3asnYtwfTFF6nz45w6lvh46jT+8AdLG8nJBOmOHckoahvr15Npn3+ephHVpTDx8fR3VJSljSNHyLRqlSbO6WSsX09XuWrxCc7p5MTH0/Ovv25p49gxOqnt2lly3RYW0sk/c4bMMGAAPX73LhkuPx9YsMDSxokTZK62bQnwnJP5Nm8m044aRak11fW98fE0QvrjHy1tnDpFpg0OploaahsJCfSlkCT6kqn3MuPj6Us0ezZ9yTinL9nGjTQqfPll2yMwxthhznmYY66quNzhS1XJyeSR+/fpil5VXBzwi1+UfP/XrgGffEJfUO3jmzbRhVNRkQVs2dnAv/715LbbttFV/aJFlscfPiTw//nPJbfdupW21bZx8yZ1gMHBlovaggK66EtNpe9S5870+P79lFu/bVvgpz+l79EPP1Cee87p4vX558kfZ85QZ9ioEY2yZs2ii9DCQvru6nR0fPZA6ujjixcTNHx96SLB358+9xkz6LgCAmhko16gBASU/FwB+rsq+FJ5zf9v735C7KzOOI7/HjQRJURjM9GYZBoXs3BETIIRkYi4s0GYYKpYpLgodGEwKXajdBEQXFoLtptgxTF/QVJsFoK0MlIRLIZSSluxpkJsZRpTukhB0RaeLp57et6ZTCb339z3npPvBy4z98zNmXPyPu993ue977ynp9i8XFx+8UUcAJ0/H9tyaipfnHjsWPzf7NsXbe5xb/z33oti5qGHou3LL+OA8dy5ONGU4uHTT6N97VrpqadyH3NzcXC3fXscHLrHe8uRI3HW55FHohBxj+LgyJHYPvv35+2Q1k64884oOtwjno4ejWJg794oMNzj+eHDcVB44EDu4913Yx2HO+6I17vHAfWxYzH2hx/ONy2an4+8cs01cd//1EeK+dtvj3GnPo4fj+Joz54Yo3sUT6+9FrH49NO5j/ffj/UCbrstCsr00dKJE1GkpTNr7rE/zs5enN8WxUhXcTlIxX23pDPu/om7fy3phKSZAfpbcPN3Kd9/uVl5SzHhq69eePP3ZGIiv0lJ+Sb0zcq72Uez8k7Wr8+Vt5RPR+7YEZV3OsJNnyk2K+8kLW6SrrZNC6Rs25Yrbym/2e3dGxv/rbdyH+vW5cpbip1LyoubnD0bz9PiJs3KO7nhhlx5SxcvbpI+Y128uMncXO7j+utz5X34cG6fno5xp0r8q680DoYel0k6rdusuKV8KrtpYmLpPtLFZikupXhj27z54tfu2LHwPuNSVPc7d8b3zYUNHnggH0ykf7NhQ5x+T2uJSxHzjz4ap/Fefz2333tvbPt0mn3Nmmh78snYL9OiPTMz8aZz332xf6b4efbZOIBeuzZ/Jr5/f/z+W26J52l/2rYt9rEk7SNLnTZPb6YXLsRZo7RwxjvvxAHw7GycbXvxxWh/7rmowF54QXrppWgbowVzhhqb110XB9ITE5Eokq1bc+UtRXI2y5V381jg2msXVt7J5GSuvKU4SDCLOEuVd5IWN0mVd9Jc3ETKX9PiJmnRj7S4yeOPRxFw8mTuY9OmXHlLeTzNyjv1sXr1wso72bgxV95SzivNyrvZR1rc5I03ch8335wrbynvv83KW4p9b9WqhZV3smFDrrylvLhJX9y9r4ekb0t6ufH8u5J+usTrvi/ptKTTk5OTvpyDB9Nx3cLH/ff31k4fw++jl74PHszbVNLpfmOMuGx/W65kH7t2dd/Hzp1Lt09Pd99Hm3HZbWz2EpfLxWaJ8XCl9tFPXK74H224+yFJh6Q49bPca5unyy59iqv7dvoYfh+99j2uiMvxHx9xuXxcSpePzXHfDvTRn0FOlX8maUvj+eZOG9Am4hLjitjEUAySuD+QNGVmt5rZakmPSTp1mX/TtUtdUdpLO30Mv49e+24BcTkGfY97Hy0ZeWyO+3agj/70fVW5JJnZbkk/kXSVpFfc/fnlXt/t1buoR0tX7xKXWFYbcdn5vV3HJnF55ek2Lgf6jNvd35T05iB9AMNGXGJcEZsYhoFuwAIAAEaLxA0AQEFI3AAAFITEDQBAQUjcAAAUhMQNAEBBSNwAABSExA0AQEFI3AAAFITEDQBAQUjcAAAUhMQNAEBBSNwAABSExA0AQEFI3AAAFITEDQBAQUjcAAAUhMQNAEBBzN1H98vMzks6u8SP1kv658gGMnq1z0+69By/6e4Tox5ML67guJTqnyNxWaba5zhQXI40cV9yEGan3f2utsexUmqfn1TnHGuc02K1z7HG+dU4p8Vqn+Og8+NUOQAABSFxAwBQkHFJ3IfaHsAKq31+Up1zrHFOi9U+xxrnV+OcFqt9jgPNbyw+4wYAAN0Zl4obAAB0gcQNAEBBWk3cZvagmX1kZmfM7Jk2xzIsZvaKmX1uZn9stN1oZr8ys487X9e1OcZBmdkWM5szsz+b2Z/M7ECnvYp5EpdlIi7LQ1z2N8/WEreZXSXpZ5K+JWla0nfMbLqt8QzRq5IeXNT2jKS33X1K0tud5yX7r6Qfuvu0pHsk7etsu+LnSVwWjbgsz6siLnueZ5sV992Szrj7J+7+taQTkmZaHM9QuPtvJP1rUfOMpNnO97OS9ox0UEPm7vPu/rvO9/+W9KGkTapjnsRloYjL8hCX/c2zzcS9SdLfGs//3mmr0U3uPt/5/h+SbmpzMMNkZlslbZf0W9UxT+KyAsRl0WrYXksaVlxycdqIefz9XRV/g2dmaySdlPQDd7/Q/FlN87wS1LS9iMt61LS9hhmXbSbuzyRtaTzf3Gmr0Tkz2yhJna+ftzyegZnZKkUQHnX3X3Saa5gncVkw4rIKNWyvBYYdl20m7g8kTZnZrWa2WtJjkk61OJ6VdErSE53vn5D0yxbHMjAzM0k/l/Shu/+48aMa5klcFoq4rEYN2+v/ViQu3b21h6Tdkv4i6a+SftTmWIY4p+OS5iX9R/E51PckfUNx1eDHkn4t6ca2xzngHHcpTuv8QdLvO4/dtcyTuCzzQVyW9yAu+5sntzwFAKAgXJwGAEBBSNwAABSExA0AQEFI3AAAFITEDQBAQUjcAAAUhMQNAEBB/ge9ncGQzYu6fAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% EMD\n",
    "G1 = ot.emd(a, b, M1)\n",
    "G2 = ot.emd(a, b, M2)\n",
    "Gp = ot.emd(a, b, Mp)\n",
    "\n",
    "# OT matrices\n",
    "pl.figure(3, figsize=(7, 3))\n",
    "\n",
    "pl.subplot(1, 3, 1)\n",
    "ot.plot.plot2D_samples_mat(xs, xt, G1, c=[.5, .5, 1])\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "# pl.legend(loc=0)\n",
    "pl.title('OT Euclidean')\n",
    "\n",
    "pl.subplot(1, 3, 2)\n",
    "ot.plot.plot2D_samples_mat(xs, xt, G2, c=[.5, .5, 1])\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "# pl.legend(loc=0)\n",
    "pl.title('OT squared Euclidean')\n",
    "\n",
    "pl.subplot(1, 3, 3)\n",
    "ot.plot.plot2D_samples_mat(xs, xt, Gp, c=[.5, .5, 1])\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "# pl.legend(loc=0)\n",
    "pl.title('OT sqrt Euclidean')\n",
    "pl.tight_layout()\n",
    "\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset 2 : Partial circle\n",
    "--------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAADSCAYAAADJ2Y62AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGNVJREFUeJzt3Xu0XGV5x/Hvr9wU0KqcFAUSooLULEurHBFX7aoHaRvQBYqaQqsLFArSUnuhZmE1yUmoVaDivVHEC4I1phQr9VJAT6itRc2hBYUiNVAQImCI5SZUiDz9Y78nGYYzM/uc2TP7Mr/PWntl9uXs/c6eyfvM++732VsRgZmZWZP8QtkFMDMzK5qDm5mZNY6Dm5mZNY6Dm5mZNY6Dm5mZNY6Dm5mZNY6Dm9kcSTpR0r+VXY4iSQpJB6TXH5W0oqD9LpL0oKSd0vxVkk4uYt9pf1+VdEJR+7PmcHCzOZP0Mkn/Luk+ST+R9E1JLy67XFUgaXEKFDuXWIZbJR0x37+PiLdExFlFHCcifhgRe0bEz+dbnpbjTUq6uG3/R0bEhf3u25qntP+AVk+Sngp8CTgNWA/sCvwG8LMBHGvniNhW9H7L1tT31W5U3qdVk1tuNlfPA4iIz0XEzyPi4Yi4IiK+CyDpFyS9U9Jtkn4s6TOSfjGte7mkO1p31vrrP/0yv0TSxZLuB06UtJOkv5R0s6QHJF0jaWHa/pclXZlajzdJWtap0JLeJOnGtI9bJJ3asu7lku6QdEYq852S3tSyfi9Jl0m6X9J3gOd2OT/fSP/em7rjXpq6Mb8p6X2StgKTkp4raUrSVkn3SPqspKe1HPNFkv4zlffvJX1e0l+1rH+VpGsl3Zta0Qen5RcBi4B/Ssdf3uF8vC29zx9JenPbuk/PHEvSmKQvpeP8RNK/ps/4CcdpabWeJOmHwFSHluxzJX0nnc8vSnpG6+fQVpZbJR0haSnwl8DvpuNdl9Zv7+bs8d2bKccJkn6Yzvk7Wo5zqKTpVKa7JZ3X5TO2OogIT55yT8BTga3AhcCRwNPb1r8Z2AQ8B9gTuBS4KK17OXBH2/a3Akek15PAo8CryX54PRl4G/A94CBAwK8CewF7ALcDbyLrgXghcA+wpEO5X0kWlAT8JvAQ8KKWcm0D1gC7AEel9U9P69eRtVL3AF4AbAb+rcNxFgMB7Nyy7MS0/z9OZX0ycADwW8BuwAKyoPj+tP2uwG3An6TyHAs8AvxVWv9C4MfAS4CdgBPSedyt/Zx2KONS4O70XvYA/i6V+YC0/tMtx3o38NFUjl3IWuma7Tgt7/0zab9Pbj8fwFXp/M0c+x+Ai+fw/bi4bf1VwMk5vnsz5fh4KtevkvU2PD+tvxp4Y3q9J3BY2f/XPPU3ueVmcxIR9wMvY0dFsSW1avZOm/w+cF5E3BIRDwJvB45T/mtQV0fEP0bEYxHxMHAy8M6IuCky10XEVuBVwK0R8amI2BYR/0lWUb6+Q7m/HBE3p338C3AFWUU941FgTUQ8GhFfAR4EDlI2EOK1wMqI+GlEXE8W2OfqRxHxoVTWhyNiU0RcGRE/i4gtwHlkQRfgMLIg+MFUnkuB77Ts6xTgYxHx7chazxeSVdSH5SzLMuBTEXF9RPyULGh08ijwLGD/VJZ/jYheN6SdTOfq4Q7rL2o59gpgWTrP/crz3Vudzv91wHVkQQ6y93mApLGIeDAivlVAeaxEDm42ZxFxY0ScGBH7kf0C3wd4f1q9D1mrY8ZtZBX13uRze9v8QuDmWbbbH3hJ6i67V9K9ZJXbM2fbqaQjJX0rda3dS9Y6G2vZZGs8/vrQQ2S/4Bek8reWq/X95fW49yVpb0nrJG1OXbAXt5RnH2BzWxBp/fv9gTPa3vvC9Hd57EP+93MuWWvoitSde2aO/bd/ht3W30bWIhzrsO1c5Pnu3dXyeuYzBjiJrMv9+5I2SnpVAeWxEjm4WV8i4vtk3VgvSIt+RFb5zlhE1iV3N/BTYPeZFenX+oL2XbbN387s17huB/4lIp7WMu0ZEae1byhpN7JW3d8Ae0fE04CvkHVR9rIllX9h23vqpFOrpn35X6dlvxIRTwXe0FKeO4F9JbWWr/X4twPvanvvu0fE53qUYcad5Hw/EfFARJwREc8Bjgb+XNIrehyn1/Hbj/0oWZdyr+9Hr/12++51FRE/iIjjgV8CzgYukbRHr7+z6nJwszlRNojjDEn7pfmFwPHATDfO54A/k/RsSXuSVeKfT62i/waeJOmVknYB3kl2zambC4CzJB2ozMGS9iIbsfk8SW+UtEuaXizp+bPsY9d0nC3ANklHAr+d5/1GNoT9UrJBILtLWkJ2jauTLcBjZNd9unkKWdfnfZL2Jbu2OONq4OfA6ZJ2lnQMcGjL+o8Db5H0knRO9kjn9Clp/d09jr+ebLDOEkm7A6s6bahs4MoBKdDel8r1WM7jdPKGlmOvAS5J57nX9+NuYLGkTvVWt+9eV5LeIGlBRDwG3JsWP9btb6zaHNxsrh4gG8jwbUk/JQtq1wNnpPWfBC4iGyDxP8D/kQ2kICLuA/6QLGBtJvul/rjRcbM4j6wyvgK4H/gE8OSIeIAsQB1H9ov9LrJf3E8Ilmnbt6b9/C/we8Blc3jPp5N1X91F1kr9VKcNI+Ih4F3AN1OXYafrYKuBF5EFjC+TBdCZfTxCNojkJLKK9g1kwfxnaf008AfAh9P72UQ2aGXGu4F3puP/xSxl/CpZN/JU+tupLu/9QOBrZIH4auBvI2JDnuN0cRHZebwLeBLZZ5Pn+/H36d+tkv5jlv12/O7lsBS4QdKDwAeA47pcM7QamBn1ZGYVJunbwEcjomNgNbMd3HIzqyBJvynpmalb8gTgYOCfyy6XWV34DiVm1XQQO3LrbgFeFxF3llsks/pwt6SZmTWOuyXNzKxxHNzMzKxxKnvNbWxsLBYvXlx2MczMrEKuueaaeyKi/eYPT1DZ4LZ48WKmp6fLLoaZmVWIpFy3v3O3pJmZNY6Dm5mZNU4hwU3SJ9PDAa/vsF6SPihpk6TvSnpREcc1a4xzzoEN2V2tJifTsg0bsuVz2cbMgOJabp8muzdbJ0eS3aPuQLJnUa0t6Lhm9dArML34xbBsGWzYwOrVad2yZdnyGb22cfAz26Gop56SPen2+g7rPgYc3zJ/E/Csbvs75JBDwqwxpqYixsYipqYCHj/fvs1qVjxxXZ5t8hzDrOaA6cgTk/JslGtH3YPbl4CXtcx/HRifZbtTgGlgetGiRYM8P2bFO/vs7YFk1aq0bGoqWz7zukNgWrUq+9+4mhUREKtZEdCyn5zbdDxGr7KZ1UQtg1vr5Jab1U6XllNfgWmWY8w5QLpVZw1RteDmbkkbDTm6Defdpdhv12a3dW7ZWU1ULbi9EvgqIOAw4Du99ufgZpXTIwD01XLKE1xydnvOq+Xolp3VxFCDG9nj3e8EHiV7cu5JwFuAt6T1Aj4C3Ax8r1eXZDi4WRX103IaRsuoj2t+udabVcDQW25FTw5uVkl1Ha3Yo2y5rgmaVUDe4OY7lJi16pIrNjkJOnyCNfecxkrOYs09p6HDJ3Zst3EjrF8PExOsWgVMTGTzGzcO/W08QY+yTU5CTG1g5dha1rCClWNriakN2Xtz/pzVUZ4IWMbklpuVolfrq6ldd93ed5VbpDZycLek2Tx1CmBNruT7vV5nNiQObmbz0PXa04gOl/f1OKsSBzezTtxKmTvf+cQqwsHNrBNfX5obny+rkLzBzaMlbfTMjBRctozVrMzurJ9GElZ6xGNZup2TbufSrEx5ImAZk1tuNii+hlQcn0sbNnK23JRtWz3j4+MxPT1ddjGsqdKz0Nbccxorx9a6tdEPn0sbIknXRMR4r+3cLWnN1C3xeOYhn+vXs4o127vVZra3Oeh2Lp38bSVycLNm6vbUal9XK063c5nn6eJmA+JuSWsud5eVz5+BFczdkjbSet4H0gbOn4GVyS03ay63Gsrnz8AK5pabjTYPGilfp8/g1FM90MQGzsHN6m+2UXnr1sGxx3rQSJk6DTYBDzSxgXO3pNVfSwtBh08QUxt8p4yqc3elzZO7JW10+BZQteKBJjYMDm5We64s66XrU7/NCuLgZrXnyrJmPNjHhsDBzerPlWW9zDbQ5Nhjs0FAeASlFcPBzeqj070Kzz3Xt9Oqk+XLt18P3f45HnccXHqpR1BaYTxa0urDoyKbzSMoLQePlrTm8ajIxvKgICuag5vVhivA5vKgICuag5vVhivABvOgICuYg5vVhyvA5vIz9qxgDm5WH64Am6t9BOVMCsDy5U4NsHkpJLhJWirpJkmbJJ05y/oTJW2RdG2aTi7iuNZw7UP/ly/fvnx7hTcxsWO5NYef4m196ju4SdoJ+AhwJLAEOF7Sklk2/XxE/FqaLuj3uDYCXMGNLo+MtT4V0XI7FNgUEbdExCPAOuCYAvZro84V3MjyyFjrVxHBbV/g9pb5O9Kydq+V9F1Jl0haONuOJJ0iaVrS9JYtWwoomtWZK7jR5ZGx1q9hDSj5J2BxRBwMXAlcONtGEXF+RIxHxPiCBQuGVDSrKldwI8wjY61PRQS3zUBrS2y/tGy7iNgaET9LsxcAhxRwXGs6V3CjyyNjrU9FBLeNwIGSni1pV+A44LLWDSQ9q2X2aODGAo5rTecKbnQ5NcD61Hdwi4htwOnA5WRBa31E3CBpjaSj02ZvlXSDpOuAtwIn9ntcayAP/bdOPHLW5shPBbDq8F3/rRs/NcDwUwGsjjz03zrwyFmbKwc3qwxXYNaJR87aXDm4WWW4ArOOPHLW5sjBzarDFZh14pGzNkcOblYdrsCsk9bUgN3TqNqJCSYfSiNnnRZgbRzcrFytw/8fShXYhg1ZBQYe+m9P5LQAy8HBzcrlisrmyqNqLQcHNyuXKyqbI4+qtTwc3KxUrqhsrjyq1vJwcLNSuaKyOfOoWsvBwc3K5YrK5sqjai0H31vSynXOOdngkYmsK3JykiywbdzoUZJm9gR57y3p4GZm9eQfRiPJN062amt/vA04Edfmxmkk1oWDm5XDFZP1y2kk1oWDm5XDFZP1yWkk1o2Dm5XCFZP1y2kk1o2Dm5XCFZP1zWkk1oWDm5XDFZP1y/lu1oVTAawcHsZtZvPgPDczGznbfyhZYznPzarLOW42IKtXl10CqwoHNxs+57iZ2YA5uNnwOcfNCjQ5CVI2wY7X7p4cbQ5uNnTOcbMiTU5CRDbBjtf+Po02BzcbOue4mdmgObjZ8DnHzQZk1aqyS2BV4eBmw+fkWxsQt/5tRiF5bpKWAh8AdgIuiIj3tK3fDfgMcAiwFfjdiLi12z6d52ZmZu2GlucmaSfgI8CRwBLgeElL2jY7CfjfiDgAeB9wdr/HtRpznpuZDVgR3ZKHApsi4paIeARYBxzTts0xwIXp9SXAK6SZgbs2cpznZmYDVkRw2xe4vWX+jrRs1m0iYhtwH7BX+44knSJpWtL0li1bCiiaVZLz3MxswCo1oCQizo+I8YgYX7BgQdnFsQFxnpuZDVoRwW0zsLBlfr+0bNZtJO0M/CLZwBIbQc5zM7NBKyK4bQQOlPRsSbsCxwGXtW1zGXBCev06YCqq+jgCGzznuZnZgPUd3NI1tNOBy4EbgfURcYOkNZKOTpt9AthL0ibgz4Ez+z2u1Zjz3MxswPw8NzMzqw0/z82qy3luZjZgDm42fM5zM7MB27nsAtgIelye22mwbK3z3MysUG652dA5z83MBs3BzYbOeW5mNmgObjZ8znMzswFzcLPhc56bmQ2Y89zMzKw2nOdm1ed8NzMbEAc3K4/z3cxsQJznZuVxvpuZDYhbblYa57uZ2aA4uFlpnO9mZoPi4Gblcb6bmQ2Ig5uVx/luZjYgznMzM7PacJ6b1Yfz3cysYA5uVj7nu5lZwZznZuVzvpuZFcwtNyud893MrGgOblY657uZWdEc3Kx8znczs4I5uFn5nO9mZgVzcLPyLV++ffDI5O4pLWBigsmHlmfrnRZgZnPk4GbV4rQAMyuAUwGsWpwWYGYFcMvNKsVpAWZWBAc3qxSnBZhZEfoKbpKeIelKST9I/z69w3Y/l3Rtmi7r55jWcE4LMLMC9NtyOxP4ekQcCHw9zc/m4Yj4tTQd3ecxrcmcFmBmBeg3uB0DXJheXwi8us/92ahrTQuYZEcKwPLlfmKAmeXWb3DbOyLuTK/vAvbusN2TJE1L+pakjgFQ0ilpu+ktW7b0WTRrBKcGmNk89EwFkPQ14JmzrHpH60xEhKROTz7dPyI2S3oOMCXpexFxc/tGEXE+cD5kDyvtWXprPqcGmNk89Gy5RcQREfGCWaYvAndLehZA+vfHHfaxOf17C3AV8MLC3oE1mlMDzGw++u2WvAw4Ib0+Afhi+waSni5pt/R6DPh14L/6PK6NCKcGmNl89Bvc3gP8lqQfAEekeSSNS7ogbfN8YFrSdcAG4D0R4eBm+Tg1wMzmoa/bb0XEVuAVsyyfBk5Or/8d+JV+jmMjrFtqgK+7mVkHvkOJVVt7agBkgS2NlnR6gJnNxsHN6sfpAWbWg58KYPXj9AAz68EtN6sdpweYWS8OblY7Tg9ogHPOeeKIV183tQI5uFn9OD2g/lqumwK+bmqFc3Cz+pktPeDYY2HdOsAjKGuh5bopK1du/7Hi66ZWFEVU8xaO4+PjMT09XXYxrC5aWnM6fIKY2uAKsw5WroSzzoIVK2DNmrJLYzUg6ZqIGO+1nVtu1gyPG0HplkAtbNgAa9dmgW3tWncrW6Ec3KwRPIKyZlpa2qzxdVMrnoObNYJHUFZUy6jIx10LPffcx7es/cR1K5ivuVkz+JpbNflzsYL5mpuNlk43WD733NlbDh5FORy+FmolcXCzZpjtBssTE/C2t/k+lCXytVAri7slrflSQFtzz2msHPN9KIfO598K5G5JM9xyGJpOA0dOPdV3k7FSOLhZo3UcRbl7h8rY1+Lmp9NjiKDzw2bNBikiKjkdcsghYda3qamIsbGIqamAlvn3vnf25VNTZZe4vtI5XM0Kn0sbGGA6csQQt9ys2TqNoty2zaP4CuTuX6ucPBGwjMktNxukVasiILJWBsRqVgRky62Ds8/e3hrbfp6mprLlM6/dcrMBI2fLrfQg1mlycLOB61QZ96rER1WnLt6pqe7rzAqUN7i5W9JGU7dnwnUaHDHquXHdErI7df964IiVJU8ELGNyy80Gyl1sc+auXKsC3C1pNj8jXYk76FvFObiZ9aNbJd7ka3K+rmYVlze4+ZqbWbtu1+Og3tfkOt1JZCZ53dfVrCnyRMAyJrfcrDR5WmZ1HWnZo/U10l2yVgu4W9JsMLoGgLK77voJzHnXm5VoKMENeD1wA/AYMN5lu6XATcAm4Mw8+3Zws0rrFgD6uV7X7/p+W2ZlB2ezHoYV3J4PHARc1Sm4ATsBNwPPAXYFrgOW9Nq3g5tVVpcA0Hfw6Hd9yzbzaplVvVvVRt5QuyV7BLeXApe3zL8deHuvfTq4WWX1O1x+gOvdMrOmq1Jwex1wQcv8G4EPd9j2FGAamF60aNEgz4/ZYPTZLdjv+tYyuGVmTVRYcAO+Blw/y3RMyzaFBLfWyS03q6VhDOjIcb3PLTNrqiq13NwtaTZj0Nfc3DKzhssb3IaRxL0ROFDSsyXtChwHXDaE45pVT69E6H7XL1++/Zl025O0Jyay5WYjRFkgnOcfS68BPgQsAO4Fro2I35G0D1lX5FFpu6OA95ONnPxkRLyr177Hx8djenp63mUzM7PmkXRNRIz32m7nfg4SEV8AvjDL8h8BR7XMfwX4Sj/HMjMzy8v3ljQzs8ZxcDMzs8bp65rbIEnaAtxWdjn6MAbcU3YhasDnKR+fp3x8nvKp83naPyIW9NqossGt7iRN57noOep8nvLxecrH5ymfUThP7pY0M7PGcXAzM7PGcXAbnPPLLkBN+Dzl4/OUj89TPo0/T77mZmZmjeOWm5mZNY6D24BIOlfS9yV9V9IXJD2t7DJVlaTXS7pB0mOSGj2Ca64kLZV0k6RNks4suzxVJemTkn4s6fqyy1JlkhZK2iDpv9L/uT8pu0yD4uA2OFcCL4iIg4H/Jnsags3ueuBY4BtlF6RKJO0EfAQ4ElgCHC9pSbmlqqxPA0vLLkQNbAPOiIglwGHAHzX1O+XgNiARcUVEbEuz3wL2K7M8VRYRN0bETWWXo4IOBTZFxC0R8QiwDjim5DJVUkR8A/hJ2eWouoi4MyL+I71+ALgR2LfcUg2Gg9twvBn4atmFsNrZF7i9Zf4OGloR2fBJWgy8EPh2uSUZjL6eCjDqJH0NeOYsq94REV9M27yDrCvgs8MsW9XkOVdmNhyS9gT+AfjTiLi/7PIMgoNbHyLiiG7rJZ0IvAp4RYx4zkWvc2Wz2gwsbJnfLy0zmzdJu5AFts9GxKVll2dQ3C05IJKWAsuBoyPiobLLY7Xkp9hboSQJ+ARwY0ScV3Z5BsnBbXA+DDwFuFLStZI+WnaBqkrSayTdAbwU+LKky8suUxWkAUmnA5eTXfhfHxE3lFuqapL0OeBq4CBJd0g6qewyVdSvA28EDk/10rWSjur1R3XkO5SYmVnjuOVmZmaN4+BmZmaN4+BmZmaN4+BmZmaN4+BmZmaN4+BmZmaN4+BmZmaN4+BmZmaN8//yexB3lMqD9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAACrCAYAAACZks5mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXuUZFdd7z/fqu7pnunpmZ4Z4pjMBIImCAkqIE8BjREuEMCwFJWHGBSNuPQCC5DXXVyjIoJXCXhRNBJMeEhAUYlBr0Iel4vySCIRJUEeISEJk4TJTCedeXZX/e4fe5+qXTXV1VXd9dpVv89avfrUPufss88537P32b/zPfvIzHAcx3EcJy9Kwy6A4ziO4zjd4w244ziO42SIN+CO4ziOkyHegDuO4zhOhngD7jiO4zgZ4g244ziO42SIN+AJks6WdEfy+8uSzu5kWcfpJ5JM0ukD3N5pcZtT8fc/Sjq/k2UdpxskXSjpg3H6wZIekFRea1kn4wZc0q2SjsSTXfy9u5fbMLOzzOzaXuY5akh6qaTPDLscg0DSUyT9q6T7JB2Q9C+SHjfscm0USddKOtp0Lfx9L7dhZs8ys8t6meeokdtNea/1HHX0S23mFzdqDzT9/ex6t9mMmX3LzLaaWaVXeY4iki6V9JaN5pP7HfNzzexTwy6EM/pI2gZcCfwq8FFgE/BU4NgQylLuQwX162b23h7n6YwovdSzJAHqYpUFM1vpdjtO78m2B96O5jBLi3DgTkl/Ienbkg5K+rtV8rlV0tPi9OZ413RQ0k3A45qWPUXSxyR9R9I3Jb0imfd4SZ+VtChpn6R3S9qUzDdJL5f0tbjMH8eLqlWZypLeJOkbkpYk3SDp1DjvhyVdF+/Ir5P0w8l6L5V0S1znm5JeLOkRwJ8CT4p30ovrONy58DAAM/uwmVXM7IiZ/bOZfQlqx/UPJO2Px+nXmjRT00L83ayxv5J0Vzz2n5Z0VjLvUknvkfQPkg4BPyZpJm7vW5LulvSnkjYn6/xG1Mq3Jf3iene6VYRFSTg+6voPJd0Wy/6ZtBzJOrXeWfOxAp7dtOx2SZfE8t8p6S2KIVFJ3yvpakn3xvU/JGkhWfdWSa+V9KVYno9Imm2zf78s6eao65skPSamPyKWeVHhUdhPJOucG5ddiuV7raQ54B+BU1TvWZ6yjkM+KDaq52sl/a6kfwEOAx8g3AC8W+uMZqqpB9+sPUlnSfqkQrTgbklvapFHc139UEn/N56rTwIPalr+iQpRiEVJ/67kkaekX0i0cYukX0nmnS3pDkmvkXRP1OovtNm3VduMqMGvx/26otCNAhfF/O+X9B+SHinpAuDFwOu00WiZmWX5B9wKPG2VeRcCH0x+nwYYMBV/fwL4CLADmAZ+NKafDdzRahvA24D/B+wETgX+s1iWcCN0A/A/CXfC3wPcAjwjzv8h4ImEiMdpwM3Aq5LtGOFuegF4MPAd4Jmr7NtvAP8BfB/hrvkHgV2xXAeBl8TtvDD+3gXMAfcD3xfzOBk4K06/FPjMsM/nAPSyDbgXuAx4FrCjaf7Lga/Ec7sTuKZJMw16a6GxXwTmgRngncCNybxLgfuAJ0etzAIXAVfEbc0Dfw/8Xlz+mcDdwCPjufvLWJbTV9m3a4FfWmXeCec3zQv447j+HqAM/HDch9Oa9r+2jQ6O1d8CfxbL/l3AF4BfifNOB54et3ES8GngnU3X3BeAU2LeNwMvX2Xffhq4k3AzrZj3QwjX9NeBNxGux3OApUT/+4CnxukdwGNaXf+j/MfG9Xwt8C3gLEJ9Md1OR3GdBk2spcNUewSN7wNeQ9D/PPCE5muphe4+C7wj6uVH4nkslt0Tj8G5hOvq6fH3SXH+s4Hvjdr4UcKNSnquV4Dfjvt+bpy/Y5V9W63NOAfYDzwmlvF/A5+O855BaBcWYhkeAZyc1Alv2bAOhi3EDQj4VuABYDH5++VmQTSLgtB4VVudKNo34LeQNKrABdQb8CcA32rK643AX6xS9lcBf5v8NuApye+PAm9YZd3/As5rkf4S4AtNaZ8lXERz8fj8FLC5aZmXMgENeNzXR8QL54548V4B7I7zriZpKID/RhcNeNN2FuK62+PvS4H3J/MFHAK+N0l7EvDNOP0+4G3JvIexdgN+uOla+J3Vzm+RF6HSOwL8YIs8a9dMso2iAV/1WAG7CWHczcn8FwLXrFL25wFfTH7fCvxc8vv3gT9dZd1/Al7ZIv2pwF1AKUn7MHBhnP4W8CvAtqb1ziaTBrwHer4W+O0WOuqkAV9s+ntEq/VpbMBfmJ7npnwvpEUDTujMrABzybJ/mSz7euADLTRx/irb+btCL/FcHyG5GQHuAZ7YYr12bcYlwO8nv7cCy3E/zgG+Sui8lZrWu5QeNOC5h9CfZ2YLyd+fd7DOqcABMzvY5bZOAW5Pft+WTD+EEHpbLP4Id/+7ASQ9TNKVCiHW+4G30hQKIlQ4BYcJQlit/N9YpXy3NaXdBuwxs0PAzxLuyvdJ+oSkh6+2o+OKmd1sZi81s72E3u0phN4ytD+/bYnhyrcpPNa4n9AIQeM5TvM+CdgC3JDo5f/E9PWW5RVN18KbO1jnQYTeUCs9tWOta2GaoLNi3/6M0BNH0m5Jl8fQ9f3AB+nPtXC7mVWbyrgnTv8Uocd1WwzPPmmV/EeaHuj59hZpnfCgJq3d3ME6q52rdpwCHIz1V0Gz1n66qd59CqHBRdKzJH0uhrYXCec81dq91vgsfzWttWszGupdM3uAEAXYY2ZXA+8mRLnukXSxgnehZ+TegK/GIUIFWfDdyfTtwE4lz906ZB/hRBY8uCnPbzaJet7Mzo3z30MIZ51hZtsIjXs3ppGU2wlhoWa+TRB0yoMJIUbM7J/M7OkEcX8FKG52bJ3lyBoz+wrhLviRMand+YX2mnoRcB7wNGA74e4bGs9xepz3E+7+z0r0st3MispjrbJ0Q0O5JaXl3g8cpbWe2rHWtXCMxkp+m5kVnoC3Eo7F98dr4efoz7VwqqS0fkuvhevM7DzCTcXfESJekPG1sA49w4n7u9H9X6ve/Z4u89sH7Ij+hIJmrX2gqd6dM7O3SZoBPgb8ASEqsQD8A+vTWrs2o6HejWXdRV1rf2RmPwScSYik/UZctCdaG9cG/EbgRxTeKdxOCGcDYGb7CGaVP5G0Q9K0pB/pIM+PAm+M6+wF/nsy7wvAkqTXK5iCytGsUBjd5gnPoB+IPd9f3cC+vRf4HUlnRJPED0jaRRDnwyS9SNKUwqsdZwJXxl7PeVFcxwiPHoreyd3AXiWmunFE0sOjYWVv/H0qIaz3ubjIR4FXSNoraQfwhqYsbgReEPXyWOD5ybx5wnG9l1CBvbVdWWLP8M+BiyQVPdM9kp6RlOWlks6UtAX4zfXtNQD/Dpwl6VEKZrALm8rxPuAdCibMsqQnxcqvHaseq3h9/TPwh5K2SSopGNd+NC4yT9DffZL2UK/Q1sN7gddK+qF4LZwu6SHA5wm9qdfF83U28FzgckmbFAyc281smXBdptfCrlhnjDQ90HMr7qb7RjblRuAnJW1RMEm+LJl3JXCypFcpGDjnJT2hXWZmdhtwPfBb8bw9hXAeCz4IPFfSM6J2ZxXMaXsJ3ocZgp9oRdKzCI8RumaNNuPDwC/E62uGcO1/3sxulfQ4SU+QNE24uTlKo9Y2cqyB/Bvwv1fj+4h/C2BmnyQYDr5EMBFc2bTeSwjPKb5CeO7xqg629VuEUMk3CRXUB4oZFl4Jeg7wqDh/P6FyKSqC1xJ6aUuEivsjXe9pnXcQLs5/JlQ+lxCeN94by/AaQkPyOuA5ZrafcJ5fTbhbPEAwdBQ3EVcDXwbukrR/A+UadZYIXoXPKzjBP0cwIr4mzv9zwvOzfwf+DfibpvXfTOjtHSRo4S+Tee8naONO4CbqlWg7Xk8wWn0uhpI/RTAmYmb/SAiFXh2XubqD/N7ddC3cEPP6KsGo8ynga0DzO/+vJZgiryNo4+2sXS+sdax+nlCB3kQ4Xn9NDGsSjt1jCKa+T7RYt2PM7K+A3yWciyVCb3qnmR0nVPTPIlyLfwL8fOylQrj+b43H/eUER3DRi/0wcEsMyY6yC32jem7Fu4DnK7is/6jNcotNWnt1TL8IOE5onC4DPlSsYGZLBJPZcwmPSL4G/FgHZXpR3M8DhBvZ9yd53k6IfL2J0FDfTrghLMXtvYJQVx6M+VzRwfZWo2WbYeE15jcTevv7CHXEC+I62wjn4SChfrgX+F9x3iXAmVFnLd+C6gTFB+qO4yRIOo1wMzZt/s6rkzmu5/Ek9x644ziO40wk3oA7juM4ToZsqAGX9ExJ/6UwCk0nJgnHGRgb0aeZ3Wpm8nCj0w8GXXe6nseTdT8DVxga8asEY8IdBBPMC83spt4Vz3HWh+vTGVVcm06v2EgP/PHA183sluj6vJzgCHScUcD16Ywqrk2nJ2zka2R7aBzJ5w6C3X9VNpVmbXN5Hsr1T73aVJi2cv39+mK6mnwR1krFvCTDUuM8AMohoqBSPbJQKoVX78pJ2pRC2lSpcmKakjROTCvH5crJu/ilOF1OvkGi2n81/E7TWmFJvtYirXiRsJpETyoxv0pyT1aJB2YlOUArlGNa/UCuVFssF9Mq1XpatRrLXE3KHqeVjHlVHKqjd92x38xOYjh0pc9NmrFZ5hpOUu0clRKBlVqlhekGDcflGnUd/ieHtJbWSsO1/0ApaneqXD/Q01G704k2N9XS6pHSaULaVHKSyk16LZHqdnVtVht0GKZTHa7EdVN9Lcfp5STtuIWqZzm5yJdbaa4SpytJ+eJuJLvdOq1SzLMkLUwvHfp2NtoEmN40Z7NbdlCZTvQUa29rUU+m2inqwqIeBCgrnv8krdBHOdVJqzSKNDshrZQsV9SJqbSLNCX1ZDFfnFh30mJepxR1ZqsYc7VVHWuptmMb1LBOrGOTi7Wob6st0tL6tFgnXbeYX0nq00L7R7+xryN99v1zogpfXrkAYLa0lSct/CTsrA9oU9kRBtk5vqM+dsTxbUGRx7bVd3Z5PuzkcjImz8p8OOArc/XDbHOxAttyvJa2dUv4wt72zUdrabtmw+h8u2bqo/SdtOkBAB40vVRL21kOabumHqilLZQOx//1L/dtibXFfKl+MmbjN+mn4/8p6ldaWasHPyrJCJArsfI9mjy6OhbnH0oqpqV4NS9W6x9uWqyEQZEOVOqjA35nZR6A/cvztbT9x8L8e4/VD+6Bo2HdxSP1/A4fDueocmi6llY6FPftUH2/px4I019566s7Ho50GDRoky08QT+OpuqXRDGt2bo2tTl8qMvm6h/ssrlwjCpb68stz4djdHy+fp6L6ePz9WO1vLVYPmkEt4VzXppfrqXNzQft7po7XEvbvSXodM9s/SNye2bC9CnT9VEfv3vqPgBOKte1vr1U6DWcvy3JOD6FXltxzOplOlwN00uJXg9Uwn7fk2jurpUwHMK+5R21tDuPhTrg20fqY6bcdSiMMnnwUP3YHlmK+luqn5eppXAcp5fqx3FTvGSnl+rHcWYplGv6geSGfCmU+ap/ffNIaxMa9TmzeYFHP/UVPHBK/Tgc3RX2//i2RDvzYZ+LehBg01yoC+c21+urbbNhemHmSC1tx0zQ1vbpJG0qpk3VdbdQDtPzpXp9Ol8K62xL0raUwrGeTe6qtsRGf1PSgM/EunA6qR8LDRY3lu3qy9Uo6tGiDg1pYfvHkvp0OTbhR5MG/KiF7R6u1o/3IQvaXqrW9Xl/rG+XKvW0+4p6d6Veny4ux/p0ub7c4vEwvXi0nrZ0NNQh/3ne73Skz4004HfSOFTf3pjWgJldDFwMsE07rXLvAdLqoZhuHAas1UBQ9fu0OmqaF0a+B1hOcqw3vT2i1VErGvNq0gVIevdAy4E1WwmzIS0uN9twAxr3Mrl7plqkHaXfHE6miz1cYfVKf0isqc9mbQLYSguPT4tD2qo/0PoITLdIS8+5mv7Xc0pLkg4G3RuKHOMZLNVveLfEa6dVQz6jZH+K3ajWG3XKxXRvr7oj1G8kV2oX4FrHcWRfsllX3Tnzievg2elXjIvjcOI+p9o5zpAHWWw4DUWNkfaLi/BJkhRn1zRoaVSgs/NaWy7dVLyZmEkr8aIxb7io40ql5EgWRRghWW2kKNcBZyh8r3UTYfSZjYx04zi9xPXpjCquTacnrLsHbmYrkn6dMFxfGXifmX25ZyVznA3g+nRGFdem0ys29AzczP6B8BGNjtBUmfLCTir3HqillZv+QxpObxdKh3ZhszSEVITTBxpKh3o4vTmUDvUiJ+GdtuH0ZLl6OD3ZyyKcXk3TBhdOT/dwVMLp3erzhPUHGk5v9WgoMRrG/30PpUMtnL6F9s/Fa+H0VLZFOL2chNX7FE5fabgAO3vENip0q03bvoVjT30cM5+4rp5YC6enx2H1x4pDD6VDUqy0xigqt9QJ2zirQX9WmOi6DKWnm0qey9fC6ekr8mpegXo4PXW2DVlao6dsx3Ecx3HWpO8u9AbKZdi50NA7KXrjE2Fsa9cTh7EztmWFguu8Va97MoxtaW5ubBs1KtMKDvTExFbrjU+CsS2Z5ca2Onmq2XEcx3EmHG/AHcdxHCdDBhpCt6kylR1zLcPlE2Fsa/eOOIydsS0nhBoGcYHWoXM3to2LsS0vbKoYuCXZvxg6nwhjW7t3xGFijW3eA3ccx3GcDBlsD7wsju+Yabjna9XbHltjW7evmEHWxrasKJXCkKktDsvwe+JpTm5sa0fnxra8sHIxZGpa/rh/E2FsW/0VM5hcY5v3wB3HcRwnQ7wBdxzHcZwMGXwIfVuZNPRdBG3c2MZYG9tGnpJqXxqrMZLhdDe2dUp7Y1teWKn40lirem0SjG1dfvwEJsLYlqeaHcdxHGfCGWgPvFpu/MZ3IPSeJ8LY1oux08GNbf2gVMLmNre2OI1kTzzNyY1t7WhtbMuMsmFzFRoV1qpeG1djW3djp8NkGNu8B+44juM4GeINuOM4juNkyGBNbCVYnhet7xvc2ObGtiFSKmFzsw1J+YTT3djWKWk4PSdUMjbNHW8IadcVNgHGtm4/fpLMHmdjm/fAHcdxHCdD1uyBS3of8BzgHjN7ZEzbCXwEOA24FfgZMzu4Vl5WhuU5WPvzfm5sq+HGtrb0Sp9WFpWtMy17v/n0xNOc3Ng2CvRKn6VSlbnNxxrSil6xG9tgUo1tnax6KfDMprQ3AFeZ2RnAVfG34wyDS3F9OqPLpbg+nT6xZgNuZp8GDjQlnwdcFqcvA57X43I5Tke4Pp1RxvXp9JP1mth2m9m+OH0XsLujtUqwMt88IP+JYes6bmxzY9u66FqfVhLL843h6bzD6W5sG2G61mdZxrbZYy3nubENJtXYtmETm5lZY6kakXSBpOslXV851PvqwnHa0U6fqTaXj7s2ncHTsT7vOzzgkjk5sN4e+N2STjazfZJOBu5ZbUEzuxi4GGDmwafaytxq4/m6sa2GG9s2Skf6TLW5deepdny+RKtebd498TSn1Y1t/euJgxvbTqBrfS48/LtsYebImhm7sW1MjG0dst4e+BXA+XH6fODj68zHcfqB69MZZVyfTk9YswGX9GHgs8D3SbpD0suAtwFPl/Q14Gnxt+MMHNenM8q4Pp1+smYI3cxeuMqsH+96ax0PyO/GthpubGtLr/RpZQgh9JRxC6evbmzr3zviaY6TZ2zrlT6nVGXHTOfPwd3YBlkb2zrER2JzHMdxnAwZ6FjoKhnTW443mL/a3x26sa2GG9v6SrUEx1cdp3/ceuJpTsN4xQzc2NYdZVXZPr22ia0VbmzL0NjWId4DdxzHcZwM8QbccRzHcTJkoCH0UqnK1i3HGoJcyx2Hd9zYVsONbT3HyrC8FdZ+bDNu4fRhfPwkzXHyjG3roawqO6Y2NpiLG9sgH2NbZ3gP3HEcx3EyZKA98HLJ2L65sStR3CNPgrFtoD1xGL6xLSOsBMtdjdM/rj3xdNqNbaNCWVW2b7AHnuLGthE3tnWI98Adx3EcJ0O8AXccx3GcDBloCH1KVXbNtg7CTYKxrW/viMNoGttyomysbKvQePY7fWwzDuH07j5+Am5sGyRlqiyUe/9FMje2wUga2zrEe+CO4ziOkyGD7YGXKuyaWfu+fVyNbUN/xQwGa2zLibJRml9uKnlx9t3YNn7Gtrwoy5jv8+uZbmwbHWNbp3gP3HEcx3EyxBtwx3Ecx8mQgZvYTtrUedB4/IxtA/z4CYyAsS0fSiVjbv5oQ/C2rrVJMLZ1+/GTei5ZGtsyo0yV+dL6PmbSLW5sg6Eb2zrEe+CO4ziOkyFr9sAlnQq8H9hNuD+42MzeJWkn8BHgNOBW4GfM7GDbjanCg6aX1lXQ8TC2jeDY6dA3Y9sg6JU+p8pVds01vqZT9P3c2AZjZ2wbAL2sO0uqsm0I3xhwY9uQjG0d0snSK8BrzOxM4InAr0k6E3gDcJWZnQFcFX87zqBxfTqjimvT6StrNuBmts/M/i1OLwE3A3uA84DL4mKXAc/rVyEdZzVcn86o4tp0+k1XJjZJpwGPBj4P7DazfXHWXYQw0Robq7KzvLEAcd7Gtu4+fgKZG9sGzEb0OV2qsHtL68c7k2Bs6z6UnuaSobFtwGy07ixhbCkNz4DnxjYYpLFtXcVrh6StwMeAV5nZ/ek8MzNWefIp6QJJ10u6fulgng5QZ/RZjz5TbR47mN83zJ086EXduXggw+GJnb7TUQ9c0jRBgB8ys7+JyXdLOtnM9kk6Gbin1bpmdjFwMcAZ37/Zdk31rk+Zn7Gtu7HTIW9j26BYrz5Tbe4+c6ftmV1cc1tubEvJ19g2KHpVdz7yBzbZrCoj8d7QMIxtQ++Jw9CMbR0VaTUkCbgEuNnM3pHMugI4P06fD3x8w6VxnC5xfTqjimvT6Ted9MCfDLwE+A9JN8a0NwFvAz4q6WXAbcDP9KeIjtMW16czqrg2nb6yZgNuZp9hlWge8OPdbKysKgulwz0f/82NbRtkAMa2ftErfW4qVdgzs3YIvWDcjG29CaWnOeVjbOsXvaw7SxhbZDTu37qL1hMGaWwbzXfEYdjGthF4ouI4juM4TrcMdCz0MsZCapTqQwlG29jW7djp9fk5GttyYlornDLddjCsVXFjW0oexrbckMQmicaLq9i/YZSokf4b20bwFTMYurFtBE694ziO4zjd4g244ziO42TIQEPowYhRaQzL9qkko2hs6/7jJ63LlI2xLSOmqfDdU/dtKI+sjW19C6WnOY2OsS03SsCMSjSEZWvx2EkwtuXy8RPohbFtXZt1HMdxHCcPBmtik5gvqakHN7je+LCNbd2/Yta+TCNvbMuIKVU5qdy7PqAb21JG0NiWGUJMU246EUVvbRKMbZmMnQ69MbZ1uynHcRzHcfJhoD1wAbMqN/bWaj248X8u3v3Y6Z2XaeSfi484ZYztpQq97v/l8lx8sD3xNKdhPRfPCxF7Z9aUCEzEc/Eux05v3u5Q2MBz8XVtwnEcx3GcPPAG3HEcx3EyZMAhdJ34kL4It06Asa03Y6e3L9NIGdsyIhgsyzQaT/oTTndjGwzf2JYfpeb6swi3ToCxrdux09NtDD2UDt0b27rN1nEcx3GcfBi4iW3qhFchIhNgbOt+7PTG/Lot09CNbRlRQmzRpqaBPgpNjr+xbfg98TSnQRvbRh+hME52MtBHfSzthgUjY2Zs63rs9Hp+WRrb1pOd4ziO4zh54A244ziO42TImgFpSbPApwmx3Cngr83sNyU9FLgc2AXcALzEzNoONFwPAzUknsiYGtu6HTs9XTdLY9sA6JU+C4Nlw+cma+F0N7a5sa17ell3FjR8bjKG0yfD2Nbt2Onpujka27rMog3HgHPM7AeBRwHPlPRE4O3ARWZ2OnAQeFlXW3ac3uD6dEYV16bTV9bsN5mZUe+wTcc/A84BXhTTLwMuBN7TyUYb7yLj/wkwtnU7dno6nb2xrU/0Wp9pb6bWG58EY9tI9sTTnPIztvWj7kyp1aMTYWzrduz0ZLksjW3rWHU1JJUl3QjcA3wS+AawaGbFcbgD2LPKuhdIul7S9d+5N7/3g53RZ736dG06/cbrTqefdNSAm1nFzB4F7AUeDzy80w2Y2cVm9lgze+xJu7r70orjdMJ69enadPqN151OP+kq2Glmi5KuAZ4ELEiaineSe4E711OAehgoSRxbY1t3Hz9J07I0tg2YXuuzCEe6sW0Uwul5G9v6UXcWTIaxrcuPnyRpeRrberS4pJMkLcTpzcDTgZuBa4Dnx8XOBz7e3aYdZ+O4Pp1RxbXp9JtO+qsnA5dJKhMa/I+a2ZWSbgIul/QW4IvAJWtlZBgVqzbeMUbc2Dbmxrb+0RN9VjGO2TIzOrEXOAnGtnx64mlOI29s61nd2Snjamzrdux0yNzY1iGduNC/BDy6RfothGc6jjM0XJ/OqOLadPrNCDztcBzHcRynWwY6fpYBK1QaQjltw+ljZ2zr7uMn6bQb2/pLFeNwdbnhULULp7uxbRTC6Xkb27qh3ePHVoybsa3rj5/A+Bjb2uA9cMdxHMfJkAH3wI2jtsKsGhKBNXriyXJ5G9u6Gzsd3Ng2KKpmLFkVqsv1xHg4JsHYlndPPM1p5I1t66LT6GUr3Ng2vsY274E7juM4ToZ4A+44juM4GTLQEHoVOGZV0sBELZw+Eca27t4RBze2DYoVxIHKNJSTEHoRTndjW0bh9M6MbTlSMQMlemrz+LEVbmxj7Ixt3gN3HMdxnAwZbA/cjENVg1Jqkgj3NZNhbOt27HRwY9tgWLEy91S20rBHRW98Ioxt49YTT3M60diWG2GkwBVm0gqr6I27sW1ijW3eA3ccx3GcDPEG3HEcx3EyZKAh9ApiyaagmgQhauH08Te29eLjJ+DGtn6wbGXuWtnelBr3yo1tmYfTW+13XhiwjIHVj1stnO7GthqTZmzzHrjjOI7jZMiAe+AlFquzUEputYve+AQY27ofO715foEb23rNspXZt7xjlbmTYGybhJ74ajmNPmbGUbPG4sfeuBvbWjMJxjbvgTuO4zhOhngD7jiO4zgZ0nGwWFIZuB6408yeI+mhwOXALuAG4CVmdrxdHhUrsVjZ0pjqDcKnAAAHP0lEQVRYhNMnwNjW/cdP0un8jG2DohfaPG5T3HlsoYOtjauxrbuPn6Q5pOQVTh8MvdBnFXHUREPFVqsT3djWjnE2tnVzaF8J3Jz8fjtwkZmdDhwEXtbLgjlOF7g2nVHG9en0hY76l5L2As8Gfhd4tSQB5wAviotcBlwIvKddPitW4kBla+uZbmwbO2PbIOiVNperZb59pPk1snaMm7Ft0l4xGwy90mcVOFydglJ6bGLF5sa2jhk3Y1unh/OdwOuon51dwKJZLXZzB7Cn1YqSLpB0vaTrHzi43GoRx9kIPdHmscUj/S+pM4n0RJ8HD1RbLeJMOGs24JKeA9xjZjesZwNmdrGZPdbMHrt1x/DuhJ3xo5fanFnY3OPSOZNOL/W5Y6f7jZ0T6SQw/GTgJySdC8wC24B3AQuSpuKd5F7gzrUyWqHMd1bm196iG9vSQsX/+RnbBkDPtLlcLXHXoW3rLMY4GNu6/fgJeDh9TXqmzyrikE03Rqhr4XQ3tnXLuBjb1jyMZvZGM9trZqcBLwCuNrMXA9cAz4+LnQ98vKclc5w1cG06o4zr0+k3G+lLvh64XNJbgC8Cl6y1woqV2b/cQQ+8wI1tmRvbhkbX2qxUSxw8tNEwer7GtskbO32odK9PK7FU3dx4+Iuq0I1tGyJnY1tXTZCZXQtcG6dvAR6/4RI4Tg9wbTqjjOvT6QfujHAcx3GcDBnox0xWqiX2H1vlPfC1cGNbWqj4f9SNbflQrZQ4sjTbwxzzMrZ1//ET8HD64KhQ4v5qkz6Lw+/Gtp6Qo7HNe+CO4ziOkyGD7YFbiXuPzW0sEze2ZWRsy4iKYGmKI/SyFw5ubDsxh5Rh9cRzo2olliqrmCzd2NZzcjG2ZVjTOo7jOI7jDbjjOI7jZMjATWwHjm5Ze8FOcWNbWqj4f5SMbfmgKkwtlVhJTlrfwulubPNwepdUKHFf86eYW+HGtp4yLGNbp3gP3HEcx3EyZKA98Eq1xOKRXvdqcGMbORjbRhtVYHpJNN4NhxPnxrbxM7blxoqVOLDShQF43IxtI9DVHKSxrVNG4LA4juM4jtMt3oA7juM4ToYMNIRerYrDh/scXh1lY1ufQukwqsa2fFAVNi1B6/c03djWKqdxMbblQMVKLC6v0wA8Fsa20XtHHPppbOsM74E7juM4ToYMtAdOVVQOTXN4ENtyY9vQjW05EUxsxlrRCje21cna2JYZoQe+wc/dZm1sG71XzKCfxrbOGIFD4DiO4zhOt3gD7jiO4zgZMvAQeulQucFeM9Bw+rCNbQMMpcPwjW05oQrMLFXp9HGDG9ta55SlsS0DVqzE4vENhtBTsjO25fHxE+iVsa0zvAfuOI7jOBkiM1t7qV5tTPoO4YZ9/8A22j8eRP770e99eIiZndTH/HtG1OZt+HkdJfq5H9loE7zuHFGGrs+BNuAAkq43s8cOdKN9YBz2Yxz2odeMwzEZh32A8dmPXjEux8P3o3d4CN1xHMdxMsQbcMdxHMfJkGE04BcPYZv9YBz2Yxz2odeMwzEZh32A8dmPXjEux8P3o0cM/Bm44ziO4zgbx0PojuM4jpMhA23AJT1T0n9J+rqkNwxy2+tF0qmSrpF0k6QvS3plTN8p6ZOSvhb/7xh2WddCUlnSFyVdGX8/VNLn4/n4iKQcPyTWE3LUJrg+J4Uc9TlO2oTR1OfAGnBJZeCPgWcBZwIvlHTmoLa/AVaA15jZmcATgV+L5X4DcJWZnQFcFX+POq8Ebk5+vx24yMxOBw4CLxtKqYZMxtoE1+fYk7E+x0mbMIL6HGQP/PHA183sFjM7DlwOnDfA7a8LM9tnZv8Wp5cIJ3APoeyXxcUuA543nBJ2hqS9wLOB98bfAs4B/jouMvL70Eey1Ca4PieELPU5LtqE0dXnIBvwPcDtye87Ylo2SDoNeDTweWC3me2Ls+4Cdg+pWJ3yTuB11AcV3gUsmtUGOc7ufPSQ7LUJrs8xJnt9Zq5NGFF9uomtQyRtBT4GvMrM7k/nWbDyj6ydX9JzgHvM7IZhl8XpD65PZ1TJWZsw2voc5NfI7gROTX7vjWkjj6RpggA/ZGZ/E5PvlnSyme2TdDJwz/BKuCZPBn5C0rnALLANeBewIGkq3kVmcz76QLbaBNfnBJCtPsdAmzDC+hxkD/w64Izo3NsEvAC4YoDbXxfxWcclwM1m9o5k1hXA+XH6fODjgy5bp5jZG81sr5mdRjjuV5vZi4FrgOfHxUZ6H/pMltoE1+eEkKU+x0GbMOL6NLOB/QHnAl8FvgH8j0FuewNlfgohxPMl4Mb4dy7hGchVwNeATwE7h13WDvfnbODKOP09wBeArwN/BcwMu3xDPC7ZaTOW2/U5AX856nPctBn3aaT06SOxOY7jOE6GuInNcRzHcTLEG3DHcRzHyRBvwB3HcRwnQ7wBdxzHcZwM8QbccRzHcTLEG3DHcRzHyRBvwB3HcRwnQ7wBdxzHcZwM+f+FjXrd9xdqlgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 50  # nb samples\n",
    "xtot = np.zeros((n + 1, 2))\n",
    "xtot[:, 0] = np.cos(\n",
    "    (np.arange(n + 1) + 1.0) * 0.9 / (n + 2) * 2 * np.pi)\n",
    "xtot[:, 1] = np.sin(\n",
    "    (np.arange(n + 1) + 1.0) * 0.9 / (n + 2) * 2 * np.pi)\n",
    "\n",
    "xs = xtot[:n, :]\n",
    "xt = xtot[1:, :]\n",
    "\n",
    "a, b = ot.unif(n), ot.unif(n)  # uniform distribution on samples\n",
    "\n",
    "# loss matrix\n",
    "M1 = ot.dist(xs, xt, metric='euclidean')\n",
    "M1 /= M1.max()\n",
    "\n",
    "# loss matrix\n",
    "M2 = ot.dist(xs, xt, metric='sqeuclidean')\n",
    "M2 /= M2.max()\n",
    "\n",
    "# loss matrix\n",
    "Mp = np.sqrt(ot.dist(xs, xt, metric='euclidean'))\n",
    "Mp /= Mp.max()\n",
    "\n",
    "\n",
    "# Data\n",
    "pl.figure(4, figsize=(7, 3))\n",
    "pl.clf()\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "pl.title('Source and traget distributions')\n",
    "\n",
    "\n",
    "# Cost matrices\n",
    "pl.figure(5, figsize=(7, 3))\n",
    "\n",
    "pl.subplot(1, 3, 1)\n",
    "pl.imshow(M1, interpolation='nearest')\n",
    "pl.title('Euclidean cost')\n",
    "\n",
    "pl.subplot(1, 3, 2)\n",
    "pl.imshow(M2, interpolation='nearest')\n",
    "pl.title('Squared Euclidean cost')\n",
    "\n",
    "pl.subplot(1, 3, 3)\n",
    "pl.imshow(Mp, interpolation='nearest')\n",
    "pl.title('Sqrt Euclidean cost')\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dataset 2 : Plot  OT Matrices\n",
    "-----------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXmYHFW5/79vIAszySQkjcoWIqvBe/UqiaI/NBmCsoiXRRyBK0kUHMCrchWCBA2ZBCVMYFx+KqAiGRYhBC4KsoQEakRwYyb3gmCGQAiELCwTlpDJJJNM5r1/nOqZmk4v1dO1nOr+fp7nPDNdXV1dy6ffc+qtU6dEVUEIIYSQZDEk7hUghBBCSPGwAieEEEISCCtwQgghJIGwAieEEEISCCtwQgghJIGwAieEEEISCCtwyxCRZhH5gfv/p0RklZ95CQkaEWkQkdsi/k76T0JHRFREDnX/v0FE5viZ1zYqqgIXkZki8oyIdInIayJyvYiMcd+7QUQ63bJDRHZ6Xj+UZVlTRaTXM0+6fCKo9VXVx1X1iKCWR/wTpCvlCP1PDra67Dq0vsA8ze56eR17Osj1UNULVPXKIJcZFRVTgYvIxQAaAcwCMBrA0QAOArBcRIa5B3Gkqo4EcBWAO9OvVfXEHIvd6JknXf4ayQaR0AjJlVgRkT1DWCz9txxbXS7Sx4UZjn04rPVKGhVRgYtIDYB5AL6pqktVdaeqvgygDsAEAF8O4TtfFpHjPK8HpCNF5BgR+YuIvCMi60RkZpZlDGihishHROR/RGSLiNwJYETG/CeLyFPuMv8iIh/yvHeZiLzofnaliJzmeW+miDwhIteKyNsi8pKIWFkRhU0YrohISkTud4/LWyLyuIgMcd8bcExFZLEnhTxTRJ7IWJY39fc5EflfEXnXdajBM98Ed95zReQVAI47/WiPd0+LyFTPZ94vIo+567IcQKrYbfUsi/7HjGUuTxWR9SLyXRF5DcAdAB4CsJ/nzHq/ItdltzN4r3cisoeIXO457itE5MAsyxlwKUZEZonIqyKyUUS+mjHvcNeTV0TkdTEZjL3c9/Z2902H69H9InKA57N/FJErReTP7vosE5FB/8aACqnAAXwS5sd+j3eiqnYCeBDAZ6JcGRE5CEbenwHYB8C/AXiqwGeGAfg9gFsBjAVwF4AveN7/CICbAJwPYByAXwK4T0SGu7O8COBTMK3weQBuE5F9PV/xcQCrYIL2QgC/EREpaUOTSRiuXAxgPcyxfi+AywFooWPqg60ApgMYA+BzAC4UkVMz5pkCYCKA40VkfwAPAPiB+32XAPhvEdnHnfd2ACtgHLgSwIziNtMf9D8ybHP5fe57B8F4eyIGZnE2DmJ98vEdAGcBOAlADYCvAujK9wEROQHmd/EZAIcBOC5jlqsBHA7j7KEA9gdwhfveEACLYLZvPIBtAH6e8fmzAXwFwHsADHO/a9BUSgWeArBJVXuyvPcqBn+msZ/bEvWWah+fOxvAI6p6h9sqflNV8wYwmNTXUAA/cT9zN4BWz/v1AH6pqn9X1V2qejOAbvdzUNW7VHWjqvaq6p0AXgDwMc/n16rqr1V1F4CbAewL8wOtNMJwZSfM/jzIPXaPq3kIQaFjmhdV/aOqPuMe03/AnNVMyZitQVW3quo2mDOuB1X1QfczywG0AThJRMYDmAxgjqp2q+qfAPyhwCrQf7uxzeVeAHNdv7YV8Z2XZDh2s8/PnQfg+6q6Sg1Pq+qbBT5TB2CRqj6rqlsBNKTfcBt09QC+rapvqeoWmMsOZwKA6/F/q2qX+94PsfvvcZGqPu9u/xKYhsCgqZQKfBOAlGS/7rKv+/5g2KiqYzLKVh+fOxDmjKAY9gOwwf2xpFnr+f8gABd7RXe/Zz8AEJHpnvTiOwD+BQN/wK+l/1HVdCt1ZJHrWA6E4co1AFYDWCYia0TkMnd6oWOaFxH5uIi0uCm7zQAuwO5BeZ3n/4MAfDHDkWNgtms/AG9n+FtoXei/3djmcoeqbh/Ed16b4ZjfzNBgPfP+ZrzbsA+AKgArPB4tdadDRKpE5JcislZE3gXwJwBjRGQPzzJe8/zfhRIdq5QK/K8wrfHTvRNFZCRMGufREL5zK8zBTvM+z//rABxS5PJeBbB/RlpvfMYyf5ghepWq3uGmLH8N4BsAxqnqGADPAkh6ijAMAndFVbeo6sWqejCAfwfwHRGZhsLHdIBDIuJ1CDAp7/sAHKiqowHcgN2PqTegrgNwa4Yj1ap6tbsue2ecQY/H4KH/8WOTy8BAF7O9LpbM38cecCtTl8F65r1O7t2GTTBp8Q96HButpgMgYC4vHAHg46paA+DT6VUrch18UxEVuKpuhrnu9TMROUFEhorIBJgUxnqY6zZB8xSAM93vmgTgDM97vwVwnIjUicieIjJORAqlUv4KoAfAt9xlno6BKcBfA7jAPSsTEakW08lpFIBqmB9LBwCIyFdgzkBIBmG4IqZz1aFucNsMYBdMOrHQMX0awAdF5N9EZAQ86TyXUQDeUtXtIvIxmNR0Pm4D8HkROV5MB58RYjoCHaCqa2HS6fNEZJiIHAPg88Vuqwf6HzOWuZyN1wGME5HRxa6Hy/MARrjHeSiA7wMY7nn/RgBXishhrhMfEpFxBZa5BMBMETlSRKoAzE2/oaq9MJ79WETeAwAisr+IHO/OMgqmgn9HRMZ6PxsWFVGBA4CqLoTpcHEtgHcB/B2mhTZNVbsHuVhvD8p0SXfcmAPT+nsb5kd0u2ddXoHpWHExgLdggl3eWyNUdQdMS3qm+5kvwdM5RVXbAHwNptPE2zBprpnueysBNMH8yF4H8K8A/jzIbS57QnDlMACPAOiEOQbXqWqLj2P6PID57mdfAPDEwMXi6wDmi8gWmI40Swps1zoAp7jb1uFu0yz0x4GzYTpzvQUTfG4psF3033JscTnHuj0H029jjZuSztUL/dIMxza5n98M8xu4EcAGmDNyb6/0H8H8JpbBbPtvAOxVYJ0eAvATmLs2Vrt/vXzXnf43N03+CMxZN9zP7QVzpv43mPR6qMjASxaEkDgRkWYA61X1+3GvCyGlQJfDp2LOwAkhhJByghU4IYQQkkCYQieEEEISSCBn4CJyk4i8ISLP5nh/qohsFnMf5lMickW2+QgJEnpJbIRekqAI6gEHzTC9P/P1Wn1cVU8uZqGpVEonTJhQwmoRW1ixYsUmVd2n8JyB0gx6SfJAL4mN+PUykApcVf/k3l8YKBMmTEBbW1vQiyUxICK+RxgLCnpJCkEviY349TLKTmyfEPP0o4dE5IO5ZhKRehFpE5G2jo6OCFePVCj0ktgIvSQFiaoC/x+Ywe8/DPMEot/nmlFVf6Wqk1R10j77RJ3ZIhUGvSQ2Qi+JLyKpwFX1XTWPsIOqPghgqJT4HFRCSoVeEhuhl8QvkVTgIvI+d+xcuGM2DwFQ6LFuhIQKvSQ2Qi+JXwLpxCYidwCYCvPouvUw4ygPBQBVvQHmQQYXikgPzGDvZ2Y8do6QwKGXxEboJQmKoHqhn1Xg/Z/D3DZBSGTQS2Ij9JIEBYdSJYQQQhIIK3BCCCEkgbACJ4QQQhIIK3BCCCEkgbACJ4QQQhIIK3BCCCEkgbACJ4QQQhIIK3BCCCEkgbACJ4QQQhIIK3BCCCEkgbACJ4QQQhIIK3BCCCEkgbACJ4QQQhIIK3BCCCEkgbACJ4QQQhIIK3BCCCEkgbACJ4QQQhJIIBW4iNwkIm+IyLM53hcR+f8islpE/iEiHw3iewnJB70kNkIvSVAEdQbeDOCEPO+fCOAwt9QDuD6g77WThQuBlhYAQEODO+38803xTmtpAU46afd5801fuDDcdS8vmkEv+6GXttAMetkPvRw8qhpIATABwLM53vslgLM8r1cB2LfQMo866ihNJI6jmkqpOo4C7uvRo1VragZOS6VUm5p2nzfX9OpqM11V5871fFdjYyybWQwA2jQg14op9NJDHi+33GemvfM7R3eNS+mWeU26a1xKX7nZTO9ZTi+DLPTSA+Plbvj1Mioh7wdwjOf1owAmFVqm9UI2NhohNIsgrljzMKdPrKzT0p/xMz2XvOn5LcbSQFlRXvZc1ahrmx3dNso4tW1USn9/kaO/Pc/RziozrbMqpc0zHG1oUG2eMXD6HfWOXned6rLZjnZVm+k9Y1Pacw29LLbQS3ca42VWEluBw6SM2gC0jR8/PsRdVCTZ5GtqUq2q2k2QRdPN63mYowroPMxRQLNOmzKluOmLpmeRN98PwxKSHiiT6OWbdxsP/361o10jU/pfH/bv5aRJ2adPnJh9+uyjHd1eY7zcNY5eFiqV7CXjZWFsq8CTnxLKlubxtPJ8tRJLbFHmEj0tqc0tTUsDZdl52b3U0R1jUtp6dtOAs+el33X0kUdU1zabFHlQXu7cO6VXfSa7l00nO9pLL+kl42XR2FaBfw7AQwAEwNEAnvSzTKuEVC1NkCCu6RT7A7AISwNlWXjZvbT/7LezKqUXfSgeL3td/7pHp/Sur5s0/M0zHd3qptt76WVFecl4OXgircAB3AHgVQA7AawHcC6ACwBc4L4vAH4B4EUAz/i5nqNxCZkjvbLsuMbSUjT19aZkLFdPPDF7Oifb9CJSULPQaNYtc7kxpYniCJTl7uXOZY7eNTm7lzedE6+XvamU/uKM7F7edA69LBsvVbN6tWi6o7PQyHg5SCI/Aw+jxCJkrtRPtpRO1J0kiukEYlkHjrjOdMIoNnj59E/MmW3zDEcfuMSkzW31Mn1m3jUypUs/26Rbq1P6lx/Sy6BLbBW433S5ZV6WQ7yMXbp8JW4hs12PGXCAbblNwe8PqL4+t9Ahw0BZOrseGZgqXzbb0Vdvz3Kd2VIvex91tGdsSts81+d3jEnprq/RyyBKrCl0PxWipV4mOV7GLl2+EomQGa20uXNVp8LR5Zim3tTPsuMs7rnoM4U1Fea2oDhamQyURZJxTL/97fL2cste9LLUEoeXqtmPqW0p6QGUUbyMXbp8JRIhs7XG3A4UtnZw8E2eTELU28ZAWSTucUoPsNI8w9HtI0Zrz8jy8nLn3im98wJHm2f0359OL+33suDZaxLdTGC8jF26fCWylJD3IOXq/Zg0IbP80DqrUjoVzm4dS/payCHCQFk8z93QP5DKjurR2ksvA4deDgLLrx8PioR6Gbt0+UrgQvpInSzHNJ0Kpy+drqr2pH6KIVcHjvr6/h9fhNekGCjzkHGsentVrz6+8rzcubd7G9rQan3xG03a20sviylxxEvr0+V+SWi8jF26fCVwIcs5/eOHzO1valIVUW1qCr3lzECZh4zjsvxy0/p/akb/vdWV5OX/ntOkvRBtO5teFlMYLwMmAfEydunylVBSQuWY/vFLnuENw/4xMlAWwHH6zkA7q1K69qIm60cxC4wsGYiXL2rS7qFVfdfJ6SXjZeQkIF7GLl2+ErSQc+eq9TfwR0m2/RHWNR4GytzQy4HQS3ppIzZ6Gbt0+UrJQsbYgkoM3hZ2iNd4GCg9ZHj50kuqy09s0h1D6WUfHi+3Vqd0xeR6ffX24CsOeumB8bIwmRmJkO4VZwWuGus1jEQQ4f5hoPSQsd+XnWiu+W6/il6q6m77p/MP5ha6ruE1uuamYFO39NID42V+svUJCOmuJVbgaSI6w0wkEba4GSgzcBztHu3eIjasWrsX0Ms+snjZ9YCjz3yyvu+2OjYsGS8jx09P9Yi9jF26fKVUIaO8ZlEOhLm/GCj7oZfFE9Y+o5f90MviidvL2KXLVwJvUVZy+scvIbXAGSj72bxZ9c4L+h+zSS994phnmS/HNJNSf8Chl4yX8ePus+WYZlLqTnRexi5dvlK0kJkpDscx1yfq63kNxw8hXgOr6EDp8fL731e97zuObhteo1u/TC9949lHU2H23/a9RutUOPSS8TI+MrzUmhrV0dF5Gbt0+UrRQmZWQPX1gbWIKoIQr4lXdKB0vexZbrxs/Wi99oyil0WR4ebaZke7htfockzru19+MNBLxsuSyNEIWo5pkcTL2KXLVwaVEmIKKDCCvL5T0YFSVdXpf1hH92h6WQr0kvHSRuLwMnbp8pVihWQnjBAI6AdeyYGSXoaA4/QNM7u1KqXv3ksv6aUFRBwvhyDpLFwItLT0vVSnBVfU/BSPYBquSF0PdVrQ0BDf6iWalhagrg5YsgSdGAnMng3U1aF5Rkv/+wsXxruOtuLxUhV44JIWXDqCXgaC66W4Xv5lymwM/Y86LJpOLwvCeBkenng5F/OB008HTjsNaPHs06Dd9FPLFyoATgCwCsBqAJdleX8mgA4AT7nlPD/L9dWiDKkTAdEB13cWTXf6xkKehcai9y1iONOxycuu4TXaM5JeBkKGlzvGpHTpZ5PjpYbkJuNlzATYKdCvl0GIuAeAFwEcDGAYgKcBHJkxz0wAPy922b5TQiF04ydZKCE9FHWgtMXLnXsbL7urRmvvo/QyDHofdbRrZGpQDz6JqWEZipuMlxYyyJjp18sgUugfA7BaVdeo6g4AiwGcEsByfdHQAMixtZi/6UIch0cxf/O3IMfW9s9QWwtcemlUq1O2ePfzFbgS8zddCDm21uZ0mxVeXvW28fLqrm9hyDR6GTQNDcCQabW4ptN4edXb1nsJxOgm42V0RBIz/dTy+QqAMwDc6Hl9DjJajjCtyVcB/APA3QAOzLO8egBtANrGjx8faiuHFEmyzsBj9/L1xU7gQ3+SLHi87KxK6aa77PVSA3aT8dJyQj4Dj0rGcQCGu/+fD8Dxs+xir+lw8IEQ8ezXWWjsey5w+tGChVJvllbgoXq5bVRKF59PL0Mlw8tHT27SbaPs9VJDdJPx0jLS+7a+vq+PQWeVu68D8jKIFPoGAAd6Xh/gTutDVd9U1W735Y0AjirpGz09KZdf1QosWQIAWHbcQpMCWrIEaG0t6StIBq3ufq6txZHTJwMLFgCzZ2PlLa39vS8nT457Lb1E7yXQ5+ZbD7fisvcvwUEza03v6PT+o5fBkuHlp55YgMc+YbWXQJwxs7UVzSeZ/UUvQya9b888E/dX1QEATu5aAixeHJyXfmr5fAXAngDWAHg/+jtkfDBjnn09/58G4G9+lp2zRclWZPwUmRpC9Gfg0Xvp7pfeVErv+aZxs2c53YySncv6x5j3M0Jb1F5qiG4W8pIxM0ZCipdBCXkSgOdhelZ+z502H8C/u/8vAPBPV9QWAB/ws1w/QvI6TvQMZgCImAJl9F6q6rM/47XvOEiKlxqSmwVT6IyZsRCml5GLW0zJJSRHELIAy8/Awyz5AiXdjBnPCG3ba+hlGnoZMzafgYdVeAZuKYPonFEpgXLbNtVbv+ro9hq6GTkZXt52Lr3Mtn/oZcSEGC+TNZRquiNGSwu2nmyGrGtBrRmyrq5uwBCBJESi6JyRJDydKn/+hRactrgO2/9rthl+dskSuhkVHi//sFcdenoq3EuAMdMGwoyXfmr5uMpuLcoSuuWTkPDZqkc5n+m4+6D3UXMr0xNfaNLeIm6xI8HT+6ijXdUV7qUqY6ZtBBwvY5cuX8knJNNA8VPMdbVKCJQ793YfF+rj2isJD3qZAWOmFYThZezS5SuZQrIjhoXwDJxe2ojj6K6xxssdYyrTS1W6aR08A2dr0hqKuLe03ANleuQ1emkBnksagOoj36tgLz37g27GTAjxMjmd2BYuBH70I/M86pPc563Ong2cfDI7YsSFZxSnRdNbgNpaM8pTa2vlPJPZ9VLr6nD5IfTSClwvpa0VVx7bgr9X1WLRiRXoJTuv2UUY8dJPLR9XGdCidBzV6mrVpiYFBj6fmh0xYsZHyxLleqbjernpcuPlT0+ll9bg9ktonlGhXrLzmp0EGC9jly5fydXbl6kgCylwbMo2ULrbnk6f+xm+k0RH99L8o+KVu5eMl5YSULyMXbp8xSskO2PYi59jU66Bkl7aC72klzYSpJexS5evsEWZICr8DHyrj3uOSQw4+UfFK3cvGS8tpdLOwPk0HYvJODZ9/RMq5Fpjbyr/dVYSE+6xeOrH5tj85suV5SXjpaUEGC+T0Qt94UIz7NySJWh4rBZz57rTTz+dz7G1Ac8zmefOBb5yS4U8k931ctMvlmDRy7X4z/90p9NLO3C9HHZ8LaZMAc69rbK8ZLy0lCDjpZ9aPq7S16JkazIR9PaqPvmkmmOUAcrxTMf18IVfGy/f+R29tJGODtWGhsrzkvHSfp55pjQvY5cuX8mWEuL1HDtJd8zILOmOGWUZKFXZA91yKtlLxkt7CcrL2KXLV9JCskdlMnAc1XnztGLOdOhlcujooJf00i5WrSo9M5SIa+ANDYA6LbgidT3mYw6uSF0PdVrQ0BD3mhEvK1cCBx0U91pER9rL79bQS9tJpeJeg+hgvEwG7e3A8OGlLSOQClxEThCRVSKyWkQuy/L+cBG5033/7yIyoagvaGkxw/8tcYeq5DOWraOjA9i0CZg4Ef2dZiwgVDddL5/4Jr1MApXmJeOlvezaBTz3HHDEEaV5WXIFLiJ7APgFgBMBHAngLBE5MmO2cwG8raqHAvgxgMaivuSaa8z40rWmNylqa83ra64pdfVJQLS3m78TJ8Kaln7obrpebji8FocfDnppOZXmJeOlvbz8MrB9e+nxMogz8I8BWK2qa1R1B4DFAE7JmOcUADe7/98NYJqIiO9vmDULWLAAaGnBY4/BtCQXLDDTiRW0twMHHgiMGhX3mgwgXDddL8c+3YLnnwe9JH6JxEvGS3tpbweGDgUOOaS05QRRge8PYJ3n9Xp3WtZ5VLUHwGYA47ItTETqRaRNRNo6OjrMxNravjTQPFzRlx5CbW0Aq09K5a23gNdeM61JywjMzXxe1t5AL0lRROIl46Wd9Paa9Pnhh5tKvBSs68Smqr9S1UmqOmmfffYBYFIMcmwt5m+6EFfgSszfdCHk2FprUmKVjjd9Xq7k87Kpi16SeGC8TB6vvAJs3RpMvAyiAt8A4EDP6wPcaVnnEZE9AYwG8KbfL2CvSrtpbwf23RcYMybuNdmNUN1Me3npKHpJiiISLxkv7aS9HdhzT+Cww0pfVhAVeCuAw0Tk/SIyDMCZAO7LmOc+ADPc/88A4Lj3uvnj/POBU08d2Kvy1FPNdBIrmzcDGzYAR2Z2wbGDcN10vXzyYnpJiiISLxkv7UPVVOCHHgoMG1b68vYsfYW0R0S+AeBhAHsAuElV/yki82FuRr8PwG8A3CoiqwG8BSNscWT23yiiDxwJD5vT55G4KTLwh0gvSQGi8jLvaxILGzYAW7YEGC/9jPYSV+HQgPZz002q111XeD6U4YhXqjrgcZUcSjV5lLOXjJf28fDDqvPnq27bln8+v15a14ktG+yUYSednaZDho1n31GQ9rLxXePllfSSWADjpZ2omz4/+GBgxIhglpmYCpydMuwjnT639Pp36KS9vHxv4+X3xtJLEj+Ml3by2mvAO+8EGy8TUYFzaEA7aW8Hxo0D3LtXKg/Xy523GS+fv5JeEgtgvLSSlStNV4QjjghumcmowFtbzcPoYcaNbXjMHahg8WLz8HoSOV1dZjjAiRMruH+M6+WIEcBxxwHzH6eXxAIYL60jnT6fMAGoqgpuucmowC+9FDjzTKCuDg1TWjBvnjv9nnuAyZNjXbVKZdUqI2Wlps8B9HkpX6rDJZNasHixO51ekjhhvLSOjg7gzTeD7y9U8m1kkTFgeMALgbrrOTxgjKxcaQZued/74l6TmHG9nHKK8VLrrofQSxI3jJdWsXKl+Rt0BZ6MM3CwZ6VNbN8OrFlT4elzl7SXC7ewJzqxB8ZLu2hvB8aPB0aODHa5iarA2bPSDp5/3gzIX6m3j3lJeznH9fKy0fSSxA/jpT28+SbwxhvhxMvEVODsWWkP7e3msaEHHBD3mliA66W4Xi47l14SC2C8tIYwR6tMzjXw1lYjYWsrFk0HUFuL5pOWYGZra//7l14a6ypWAjt2AKtXAx/9KNPnAAZ4efXxwIqaWrz22SU4j16SOGG8tIb2dmD//YHRo4NfdnLOwC+91HTAmDwZMx80Lcmv3GJeo66OvSsj4oUXgJ4eps/78Hg5q7UOE15qwddup5ckZhgvreCdd4CNG8OLl8k5A0/D3pWx0t4OVFebDhnEQ20t9M4lOOOUOqyll8QWGC9jJeyHPSXnDNyFvSvjY+dO04HtiCOAIYkzJ1waGoA9P1OLpi56SeyB8TJe2tvNrbZjx4az/MSFYfaujI8XXzSVeEUP3pKDtJffH8fe6MQeGC/jY8sWYN26cC83Jq4C7+tdefrpaIFJD2092e1d2dLCoQJDpL3dPEVnwoS418RCXC+HfMF4ec+X6CWxAMbL2Ag7fQ4ksQJP964880zcX1UHADi5yx3nl50zQmPXLjN86gc+AOyxR9xrYyEeL/8wog7d3fSSWADjZWy0twOpVLgPe0peJzbPrQ/V95vOGbW40Izzy84ZobFmDdDdzd7nOfF6+cASfPHzpjOb3nMPh1Yl8cF4GQtbtwJr1wLHHBPu95R0Bi4iY0VkuYi84P7dO8d8u0TkKbfcV8p3pmHnjGhpbweGDTMPo08CcbnZ0AAMmVaLa7s4tCrZnTi9ZLyMjueei+hhT6o66AJgIYDL3P8vA9CYY77OwSz/qKOO0rw4jmoqpfMwRzWVMq9J4OzapdrYqHr33YNfBoA2LcG1YkuYbvrxstf1smtkSnuW00tbqTQvGS+j4dZbVX/6U9Xe3sF93q+XpV4DPwXAze7/NwM4tcTl+cczVGAnRgKzZwN1dWie0dL/PjtoBMLatcC2bYlLn8fjpmdo1U6MxJ8+ORu9X6SXpI9YvWS8DJ9t24CXXormYU+lVuDvVdVX3f9fA/DeHPONEJE2EfmbiOQVVkTq3XnbOjo6cs+Y7pxRW4sjp08GFiwAZs/Gylta+2VlB41AWLkS2HNP4NBD416TogjUzcF4OfGcyZjylwV4/Bh6SfqI3UvGy3BZtco87CmS220LnaIDeATAs1nKKQDeyZj37RzL2N/9ezCAlwEc4ic9UDAl5IXpoVDo7VW99lrVO+8sbTkIIVUZl5vFeLnxt452VtFLW6lULxkvw+P221V/9KPBp89VA0yhq+pxqvoj9NXaAAAMbUlEQVQvWcq9AF4XkX0BwP37Ro5lbHD/rgHwRwAfKfS9xcAOGuGxbh3Q2Wln+tx2NxsagP3+g6OzVRpJ8JLxMhy6u82AV1Gkz4HSU+j3AZjh/j8DwL2ZM4jI3iIy3P0/BeD/AVhZ4vcOgKMNhcfKlea+78MPj3tNiiZ2NzOfFX5J9fXoXkovKxxrvGS8DJ4XXjBjZkQ1WmWpFfjVAD4jIi8AOM59DRGZJCI3uvNMBNAmIk8DaAFwtaoGWoGzg0Y4qJrbxw45BBg+PO61KZr43czo0Pb4MbOh9LLSscZLxsvgaW8HRo4EDjwwmu8raSAXVX0TwLQs09sAnOf+/xcA/1rK9xRkQAcNAAvqTAeNi1uBmeiTlRTHxo3Au+8mc6wHK9zM8LL23jq0HO12HJoJelmB2Ogl42Uw7NxpzsA//OFo0udAEodSzUb62bcAZt7sPj5vwQKMRGe/jEmsheJg4ULTAodJnw8ZAkx8jS3yQZHh5R53L8Gn/my83HUGvSwKj5d98ExxcDBeBovr5urVphKfOBGRuVkeFbgHdtAokcmTgbo6qNOC9nbg6G0tGH4ObzEplfTjRq/darz84Vv0sihcL/sqcd76FAiMlwHgutmxpAV77QVMeClCN/10VY+rFHVbhBfvLRJVVapNTaqqOneu5/3GxsEtuxJwHN01NqV//PQc3TEmmFtMEPGIV2GWILzsHlqlfz6jSbu76aVv3BHuNp43R3sDuvWJXirjZQD0LHd0a3VKn6sL5rY8v17GLl2+MighXRnVcczWNTWpiqg2NZnXnvdJbl6ePkcV0O7vzglkeRUfKDO8fOOyJu2F6JNn0Uu/7Nyp+o/TjJfvXEQvMwvjZXy88ILqHz9t3NQ5pbtZuRV4Y2OfbH0tyKYm1epqDlqQjRz7q7e6WjdfFNz+qvhAmWU/v3xRk3YPNV4GdUZZNmTsr507VdvONvtrw3n0MlthvIyAbPvLcbT3a/XaMzalu77HM/DShMxg7lyzlfNgWkfzMEcBz86vdPK0wDPfL4WKD5QZ0MsCZHj55FkmY7HuO/QyV6GXEZAZLx1HtaZGdfTofhcDcJMVuBde48lPjv0z4P0S9w8DZRYyron/7UtNum0bvezD81S37qFVuv479DJfYbyMiMxhaOvrd6+sS9xHrMDT8BpPXqJqcTNQZpDh5cZZ5gzz8dPppSq9HExhvAwf27yMXbp8JRAheY1nIDHtDwbKDLIchzdm918T7xlb2V6uXq36yOf69we9ZLyMBcvjZezS5SuBpYQ8ZGtBzUKjLpq+e8eEskwTxdTCZqDMTy4vbzqnsrzsfdR4ufSzJiPReSW99FsYL0PA8ngZu3T5ShhCquru1zCamnbvmFDOLc0YrnExUPrAc1y216R06WebdNuolHY9UBlevnuvuZd2HubozmFVurORXhZTGC9DwuJ4Gbt0+UooQmbrReiRsqzSRFnSP4umOzoLjZH3MmWgLECGl72POrpjTEofPqGp73niZXOrWYaXvb2q15xEL0stjJclksB4Gbt0+UooQvo8SGWRJrLox8dAWYAK9vKBWY52VqX0ybOadNc4ejnYQi9LJIHxMnbp8pXQUkLZKNc0kSXbxUA5SDy3Um2tTunS45t0e03y0+pdDzjaNbJ/u9Zf3NSXYaCXyfDShrgSOJZsFyvwYvDb8vLc72ddK9PyljID5SDI8HL7Q452u9fG+9Lq45Ll5ebNqj+YRi/DKIyXRVAm8TJ26fKVyIT0eTCnwqT6Ym1l5hjKT+vrrUn/ZIOBchAU4eWWvSw4+8lzy82rt5t1u3mm+Q0985Xo0+XZoJeDIEnxMsf6pr1MeryMXbp8JdKUUDYy0ymOs/u0qFuZuVq/2dbNorQWA2WAeI7z9pqU3nauo80zHO2q9nR2i+PsJ8PNbQ+ajnh//VJ/xmB7TUq7fkAvwyg2eWlNvPSsV8HKOoHxsiRhAHwRwD8B9AKYlGe+EwCsArAawGV+lx+rkFkOemdVSqfC8dfKDELUAmfb3h+K7fdrRh0ow3TTOi/3yu7luyNSuu3BCL1sbNSt9zvaPdq42VmV0os+5NDLSvXSlniZnu6JmYuml4eXpco4EcARAP6YS0YAewB4EcDBAIYBeBrAkX6WH6uQRVSe+VqeA0TNdQ/hiSf6Totn+1H03daQbT0sIYZAGZqbtnrZN2746JTe9XVzVu49833xRkff+Z2TvbNYtgBaX29Klu/yevnijY5ur0npZR/fPSgCagajoZcV6+Wg42WJXqaXka2yBkzaP+leBiVlPhk/AeBhz+vZAGb7WW7sKaFMfLYyc1aouVI0+VI3ftJS+X4AlkgZV6oyDDeT7OWUKaq/Pa9/wJTu0Sltv87Rlxc52jM2pRt/a5bx9j2O7hpVo7tGjdZ1t5hpz11vzq4fa3D0gUv6GwedVSaNf+utqk/92KTNbUxLZoNehkip8TJXXPM8/ctXvNQsyy4TL6OQ8QwAN3penwPg536Wa52QJaa087b6soiXbRmhpqBCxNJAOSg3k+hlbyqlb97taH19di+nTNEBZ+ydVSltnuHsNu3mmY4ef3x+t31lneglvfQRLwtV7MXE3HL00o9ojwB4Nks5xTNPYDICqAfQBqBt/Pjx4e+pUsl35uvjuks6ePqW1/KKOhdhBMoo3SxXL3tT5p7yiy/O7mBRQTXfNUhLoZcRU0S8zFcp00ufFbivhVRKCj0bQd3WlcC0eDFYeqZTHqnKbBTrZRldrikGehkx9NIXNlXgewJYA+D96O+Q8UE/y02EkLkYxL2HSUuLF4OlgXJQbpadl7kCaK5rjfSSXgYNvRxAJBU4gNMArAfQDeD1dKsRwH4AHvTMdxKA52F6Vn7P7/ITLWQ2ckmaqxd6AsXLRdSBMkw3y85L1exu5urtSy/pZVTQy7xFzLx2MmnSJG1ra4t7NUgAiMgKVZ0U93oEAb0sH+glsRG/Xg6JYmUIIYQQEiyswAkhhJAEwgqcEEIISSCswAkhhJAEwgqcEEIISSCswAkhhJAEwgqcEEIISSCswAkhhJAEwgqcEEIISSCswAkhhJAEwgqcEEIISSCswAkhhJAEwgqcEEIISSCswAkhhJAEwgqcEEIISSCswAkhhJAEwgqcEEIISSCswAkhhJAEUlIFLiJfFJF/ikiviEzKM9/LIvKMiDwlIm2lfCchhaCXxFboJgmSPUv8/LMATgfwSx/z1qrqphK/jxA/0EtiK3STBEZJFbiqtgOAiASzNoQEAL0ktkI3SZBEdQ1cASwTkRUiUp9vRhGpF5E2EWnr6OiIaPVIhUIvia34cpNeVjYFz8BF5BEA78vy1vdU9V6f33OMqm4QkfcAWC4iz6nqn7LNqKq/AvArAJg0aZL6XD6pMOglsZUo3aSXlU3BClxVjyv1S1R1g/v3DRH5HYCPAcgaKL2sWLFik4iszZicAlBO14XKbXuA7Nt0UJBfQC9Dp9y2B4jASyA+NyvES6D8tmnQXpbaia0gIlINYIiqbnH//yyA+X4+q6r7ZFlem6rm7L2ZNMpte4BkbBO9zE+5bQ+QnG0arJuV4CVQfttUyvaUehvZaSKyHsAnADwgIg+70/cTkQfd2d4L4AkReRrAkwAeUNWlpXwvIfmgl8RW6CYJElFN1mUTtr7spxy3qRDlts3ltj1AeW5TIcpxm8ttm2I7A4+JX8W9AgFTbtsDlOc2FaLctrnctgcoz20qRDluc7lt06C3J3Fn4IQQQghJ5hk4IYQQUvGwAieEEEISSCIrcL8PBLAdETlBRFaJyGoRuSzu9SkVEblJRN4QkWfjXpc4oJd2Qi/ppY0E4WUiK3D0PxCg4KAbtiIiewD4BYATARwJ4CwROTLetSqZZgAnxL0SMUIv7aQZ9JJe2kczSvQykRW4qrar6qq416NEPgZgtaquUdUdABYDOCXmdSoJd6jHt+Jej7igl3ZCL+mljQThZSIr8DJhfwDrPK/Xu9MIiRN6SWyEXmYh9KFUB0tADwQgJFDoJbERelmZWFuBB/FAAMvZAOBAz+sD3GnEYuglsRF6WZkwhR4frQAOE5H3i8gwAGcCuC/mdSKEXhIboZdZSGQFnuuBAElCVXsAfAPAwwDaASxR1X/Gu1alISJ3APgrgCNEZL2InBv3OkUJvbQTekkvbSQILzmUKiGEEJJAEnkGTgghhFQ6rMAJIYSQBMIKnBBCCEkgrMAJIYSQBMIKnBBCCEkgrMAJIYSQBMIKnBBCCEkg/wfPnU7C6ZZu4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% EMD\n",
    "G1 = ot.emd(a, b, M1)\n",
    "G2 = ot.emd(a, b, M2)\n",
    "Gp = ot.emd(a, b, Mp)\n",
    "\n",
    "# OT matrices\n",
    "pl.figure(6, figsize=(7, 3))\n",
    "\n",
    "pl.subplot(1, 3, 1)\n",
    "ot.plot.plot2D_samples_mat(xs, xt, G1, c=[.5, .5, 1])\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "# pl.legend(loc=0)\n",
    "pl.title('OT Euclidean')\n",
    "\n",
    "pl.subplot(1, 3, 2)\n",
    "ot.plot.plot2D_samples_mat(xs, xt, G2, c=[.5, .5, 1])\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "# pl.legend(loc=0)\n",
    "pl.title('OT squared Euclidean')\n",
    "\n",
    "pl.subplot(1, 3, 3)\n",
    "ot.plot.plot2D_samples_mat(xs, xt, Gp, c=[.5, .5, 1])\n",
    "pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples')\n",
    "pl.axis('equal')\n",
    "# pl.legend(loc=0)\n",
    "pl.title('OT sqrt Euclidean')\n",
    "pl.tight_layout()\n",
    "\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
